Subtracting very small probabilities - How to compute?
up vote
3
down vote
favorite
This question is an extension of a related question about adding small probabilities. Suppose you have log-probabilities $ell_1 geqslant ell_2$, where the corresponding probabilities $exp(ell_1)$ and $exp(ell_2)$ are too small to be distinguished from zero in the initial computational facility being used (e.g., base R
). We want to find the log-diffference of these probabilities, which we denote by:
$$ell_- equiv ln big( exp(ell_1) - exp(ell_2) big)$$
Questions: How can you effectively compute this log-difference? Can this be done in the base R
? If not, what is the simplest way to do it with package extensions?
r computational-statistics underflow
add a comment |
up vote
3
down vote
favorite
This question is an extension of a related question about adding small probabilities. Suppose you have log-probabilities $ell_1 geqslant ell_2$, where the corresponding probabilities $exp(ell_1)$ and $exp(ell_2)$ are too small to be distinguished from zero in the initial computational facility being used (e.g., base R
). We want to find the log-diffference of these probabilities, which we denote by:
$$ell_- equiv ln big( exp(ell_1) - exp(ell_2) big)$$
Questions: How can you effectively compute this log-difference? Can this be done in the base R
? If not, what is the simplest way to do it with package extensions?
r computational-statistics underflow
1
Possible duplicate of Computation of likelihood when $n$ is very large, so likelihood gets very small?
– Xi'an
1 hour ago
"where the corresponding probabilities exp(ℓ1) and exp(ℓ2) are too small to be distinguished from zero" - I think you mean one instead of zero?
– Don Hatch
3 mins ago
add a comment |
up vote
3
down vote
favorite
up vote
3
down vote
favorite
This question is an extension of a related question about adding small probabilities. Suppose you have log-probabilities $ell_1 geqslant ell_2$, where the corresponding probabilities $exp(ell_1)$ and $exp(ell_2)$ are too small to be distinguished from zero in the initial computational facility being used (e.g., base R
). We want to find the log-diffference of these probabilities, which we denote by:
$$ell_- equiv ln big( exp(ell_1) - exp(ell_2) big)$$
Questions: How can you effectively compute this log-difference? Can this be done in the base R
? If not, what is the simplest way to do it with package extensions?
r computational-statistics underflow
This question is an extension of a related question about adding small probabilities. Suppose you have log-probabilities $ell_1 geqslant ell_2$, where the corresponding probabilities $exp(ell_1)$ and $exp(ell_2)$ are too small to be distinguished from zero in the initial computational facility being used (e.g., base R
). We want to find the log-diffference of these probabilities, which we denote by:
$$ell_- equiv ln big( exp(ell_1) - exp(ell_2) big)$$
Questions: How can you effectively compute this log-difference? Can this be done in the base R
? If not, what is the simplest way to do it with package extensions?
r computational-statistics underflow
r computational-statistics underflow
asked 3 hours ago
Ben
21k22499
21k22499
1
Possible duplicate of Computation of likelihood when $n$ is very large, so likelihood gets very small?
– Xi'an
1 hour ago
"where the corresponding probabilities exp(ℓ1) and exp(ℓ2) are too small to be distinguished from zero" - I think you mean one instead of zero?
– Don Hatch
3 mins ago
add a comment |
1
Possible duplicate of Computation of likelihood when $n$ is very large, so likelihood gets very small?
– Xi'an
1 hour ago
"where the corresponding probabilities exp(ℓ1) and exp(ℓ2) are too small to be distinguished from zero" - I think you mean one instead of zero?
– Don Hatch
3 mins ago
1
1
Possible duplicate of Computation of likelihood when $n$ is very large, so likelihood gets very small?
– Xi'an
1 hour ago
Possible duplicate of Computation of likelihood when $n$ is very large, so likelihood gets very small?
– Xi'an
1 hour ago
"where the corresponding probabilities exp(ℓ1) and exp(ℓ2) are too small to be distinguished from zero" - I think you mean one instead of zero?
– Don Hatch
3 mins ago
"where the corresponding probabilities exp(ℓ1) and exp(ℓ2) are too small to be distinguished from zero" - I think you mean one instead of zero?
– Don Hatch
3 mins ago
add a comment |
1 Answer
1
active
oldest
votes
up vote
2
down vote
To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:
$$begin{equation} begin{aligned}
exp(ell_1) - exp(ell_2)
&= exp(ell_1) (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
This result converts the difference to a product, which allows us to present the log-difference as:
$$begin{equation} begin{aligned}
ell_-
&= ln big( exp(ell_1) - exp(ell_2) big) \[6pt]
&= ln big( exp(ell_1) (1 - exp(-(ell_1 - ell_2))) big) \[6pt]
&= ell_1 + ln (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
In the case where $ell_1 = ell_2$ we obtain the expression $ell_+ = ell_1 + ln 0 = -infty$. Using the Maclaurin series expansion for $ln(1-x)$ we obtain the formula:
$$begin{equation} begin{aligned}
ell_-
&= ell_1 - sum_{k=1}^infty frac{exp(-k(ell_1 - ell_2))}{k} quad quad quad text{for } ell_1 neq ell_2. \[6pt]
end{aligned} end{equation}$$
Since $exp(-(ell_1 - ell_2)) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $ell_1 - ell_2$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.
Implementation in base R: It is possible to compute this log-difference accurately in base R
using the log1p
function. This is a primitive function in the base package that computes the value of $ln(1+x)$ for an argument $x$ (with accurate computation even for $x ll 1$). This primitive function can be used to give a simple function for the log-difference:
logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }
Implementation with VGAM
package: Machler (2012) analyses accuracy issues in evaluating the function $ln(1-exp(-|x|))$, and suggests that use of the base R
functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp
function in the VGAM
package. This gives you the an alternative function for the log-difference:
logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "65"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstats.stackexchange.com%2fquestions%2f383523%2fsubtracting-very-small-probabilities-how-to-compute%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
2
down vote
To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:
$$begin{equation} begin{aligned}
exp(ell_1) - exp(ell_2)
&= exp(ell_1) (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
This result converts the difference to a product, which allows us to present the log-difference as:
$$begin{equation} begin{aligned}
ell_-
&= ln big( exp(ell_1) - exp(ell_2) big) \[6pt]
&= ln big( exp(ell_1) (1 - exp(-(ell_1 - ell_2))) big) \[6pt]
&= ell_1 + ln (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
In the case where $ell_1 = ell_2$ we obtain the expression $ell_+ = ell_1 + ln 0 = -infty$. Using the Maclaurin series expansion for $ln(1-x)$ we obtain the formula:
$$begin{equation} begin{aligned}
ell_-
&= ell_1 - sum_{k=1}^infty frac{exp(-k(ell_1 - ell_2))}{k} quad quad quad text{for } ell_1 neq ell_2. \[6pt]
end{aligned} end{equation}$$
Since $exp(-(ell_1 - ell_2)) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $ell_1 - ell_2$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.
Implementation in base R: It is possible to compute this log-difference accurately in base R
using the log1p
function. This is a primitive function in the base package that computes the value of $ln(1+x)$ for an argument $x$ (with accurate computation even for $x ll 1$). This primitive function can be used to give a simple function for the log-difference:
logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }
Implementation with VGAM
package: Machler (2012) analyses accuracy issues in evaluating the function $ln(1-exp(-|x|))$, and suggests that use of the base R
functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp
function in the VGAM
package. This gives you the an alternative function for the log-difference:
logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }
add a comment |
up vote
2
down vote
To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:
$$begin{equation} begin{aligned}
exp(ell_1) - exp(ell_2)
&= exp(ell_1) (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
This result converts the difference to a product, which allows us to present the log-difference as:
$$begin{equation} begin{aligned}
ell_-
&= ln big( exp(ell_1) - exp(ell_2) big) \[6pt]
&= ln big( exp(ell_1) (1 - exp(-(ell_1 - ell_2))) big) \[6pt]
&= ell_1 + ln (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
In the case where $ell_1 = ell_2$ we obtain the expression $ell_+ = ell_1 + ln 0 = -infty$. Using the Maclaurin series expansion for $ln(1-x)$ we obtain the formula:
$$begin{equation} begin{aligned}
ell_-
&= ell_1 - sum_{k=1}^infty frac{exp(-k(ell_1 - ell_2))}{k} quad quad quad text{for } ell_1 neq ell_2. \[6pt]
end{aligned} end{equation}$$
Since $exp(-(ell_1 - ell_2)) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $ell_1 - ell_2$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.
Implementation in base R: It is possible to compute this log-difference accurately in base R
using the log1p
function. This is a primitive function in the base package that computes the value of $ln(1+x)$ for an argument $x$ (with accurate computation even for $x ll 1$). This primitive function can be used to give a simple function for the log-difference:
logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }
Implementation with VGAM
package: Machler (2012) analyses accuracy issues in evaluating the function $ln(1-exp(-|x|))$, and suggests that use of the base R
functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp
function in the VGAM
package. This gives you the an alternative function for the log-difference:
logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }
add a comment |
up vote
2
down vote
up vote
2
down vote
To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:
$$begin{equation} begin{aligned}
exp(ell_1) - exp(ell_2)
&= exp(ell_1) (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
This result converts the difference to a product, which allows us to present the log-difference as:
$$begin{equation} begin{aligned}
ell_-
&= ln big( exp(ell_1) - exp(ell_2) big) \[6pt]
&= ln big( exp(ell_1) (1 - exp(-(ell_1 - ell_2))) big) \[6pt]
&= ell_1 + ln (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
In the case where $ell_1 = ell_2$ we obtain the expression $ell_+ = ell_1 + ln 0 = -infty$. Using the Maclaurin series expansion for $ln(1-x)$ we obtain the formula:
$$begin{equation} begin{aligned}
ell_-
&= ell_1 - sum_{k=1}^infty frac{exp(-k(ell_1 - ell_2))}{k} quad quad quad text{for } ell_1 neq ell_2. \[6pt]
end{aligned} end{equation}$$
Since $exp(-(ell_1 - ell_2)) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $ell_1 - ell_2$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.
Implementation in base R: It is possible to compute this log-difference accurately in base R
using the log1p
function. This is a primitive function in the base package that computes the value of $ln(1+x)$ for an argument $x$ (with accurate computation even for $x ll 1$). This primitive function can be used to give a simple function for the log-difference:
logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }
Implementation with VGAM
package: Machler (2012) analyses accuracy issues in evaluating the function $ln(1-exp(-|x|))$, and suggests that use of the base R
functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp
function in the VGAM
package. This gives you the an alternative function for the log-difference:
logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }
To see how to deal with differences of this kind, we first note a useful mathematical result concerning differences of exponentials:
$$begin{equation} begin{aligned}
exp(ell_1) - exp(ell_2)
&= exp(ell_1) (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
This result converts the difference to a product, which allows us to present the log-difference as:
$$begin{equation} begin{aligned}
ell_-
&= ln big( exp(ell_1) - exp(ell_2) big) \[6pt]
&= ln big( exp(ell_1) (1 - exp(-(ell_1 - ell_2))) big) \[6pt]
&= ell_1 + ln (1 - exp(-(ell_1 - ell_2))). \[6pt]
end{aligned} end{equation}$$
In the case where $ell_1 = ell_2$ we obtain the expression $ell_+ = ell_1 + ln 0 = -infty$. Using the Maclaurin series expansion for $ln(1-x)$ we obtain the formula:
$$begin{equation} begin{aligned}
ell_-
&= ell_1 - sum_{k=1}^infty frac{exp(-k(ell_1 - ell_2))}{k} quad quad quad text{for } ell_1 neq ell_2. \[6pt]
end{aligned} end{equation}$$
Since $exp(-(ell_1 - ell_2)) < 1$ the terms in this expansion diminish rapidly (faster than exponential decay). If $ell_1 - ell_2$ is large then the terms diminish particularly rapid. In any case, this expression allows us to compute the log-sum to any desired level of accuracy by truncating the infinite sum to a desired number of terms.
Implementation in base R: It is possible to compute this log-difference accurately in base R
using the log1p
function. This is a primitive function in the base package that computes the value of $ln(1+x)$ for an argument $x$ (with accurate computation even for $x ll 1$). This primitive function can be used to give a simple function for the log-difference:
logdiff <- function(l1, l2) { l1 + log1p(-exp(-(l1-l2))); }
Implementation with VGAM
package: Machler (2012) analyses accuracy issues in evaluating the function $ln(1-exp(-|x|))$, and suggests that use of the base R
functions may involve a loss of accuracy. It is possible to compute this log-difference more accurately in using the log1mexp
function in the VGAM
package. This gives you the an alternative function for the log-difference:
logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }
answered 3 hours ago
Ben
21k22499
21k22499
add a comment |
add a comment |
Thanks for contributing an answer to Cross Validated!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Some of your past answers have not been well-received, and you're in danger of being blocked from answering.
Please pay close attention to the following guidance:
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstats.stackexchange.com%2fquestions%2f383523%2fsubtracting-very-small-probabilities-how-to-compute%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
1
Possible duplicate of Computation of likelihood when $n$ is very large, so likelihood gets very small?
– Xi'an
1 hour ago
"where the corresponding probabilities exp(ℓ1) and exp(ℓ2) are too small to be distinguished from zero" - I think you mean one instead of zero?
– Don Hatch
3 mins ago