Generate random numbers from a particular function
in [Matlab]
|
Prev: Setting alphadata on superimposed transparent image affects underlying image.
Next: Locating GUIDE Gui Figure Handle/Passing Data
From: Gonzalo on 1 Jun 2010 01:37 "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hu22nj$jv4$1(a)fred.mathworks.com>... > "Gonzalo " <glpita(a)gmail.com> wrote in message <hu1o0g$44h$1(a)fred.mathworks.com>... > > ...... > > I got the integral of f(x) and normalized the function so that Int(f(x)) = 1; the 2nd condition is also met thanks to the type of function chosen. Now, could please explain a little bit more on how to use fzero to generate the random numbers from f(x)?? > > ...... > > Gonzalo, you still haven't stated definitely that f(x) is to be the probability density function for your generated random variable. Is that true? > > On the assumption that it is indeed true, I've looked a little more closely into the problem. To make f(x) a valid density function, one condition is certainly true and I strongly suspect another condition must also hold as well. Certainly you must have p(1) = 1/(2*(p3-p2)) in order to have an integral from minus infinity to plus infinity of one. That would be your normalization factor. The other strongly suspected condition is that p4 and p5 must be exactly equal. It can be shown rigorously that if you set these equal and also have p3>p2, then f(x) is always positive. Otherwise I am almost sure (99%) that somewhere along the x range a value can be found that makes f(x) negative, and I don't think you want that. If you lived up to both the above conditions, it would unfortunately leave you with only three independent parameters. > > However there is a decided advantage for you in doing so. Applying this constraint would leave you with a cdf function whose inverse is directly computable using the inverse hyperbolic tangent function atanh. I can show you the details later. No call on fzero would be necessary. You could then make direct use of the rand function using this inverse cdf function. > > Does this possibility interest you, given that I am only fairly certain about the above? > > Roger Stafford Mr Stafford, Yes, I need f(x) to be the pdf for my generated random variable. Now, maybe I'm confused here, but I don't have that freedom to set p4=p5. Those coefficients come from a nonlinear regression of data, and f(x) fits the data very well. Regarding p1, it is approx equal to = 1/(2*(p3-p2)) * 100. Here are the coeffs: p = [4.4782 16.8692 28.0382 6.6695 3.8967]; I've plotted f(x) and seem to be always>0. Here's the code syms x positive f = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))); % x is defined from 0 to 80. norm_c = double(int(f,-1000,1000)); % Normalization coefficient f = f / norm_c; F = int(f,x); % integral of f = 1 ord = subs(F,0); F = F - ord; Given what I wrote above, do I still need the fzero function? I'm stuck at it..
From: Roger Stafford on 1 Jun 2010 05:27 "Gonzalo " <glpita(a)gmail.com> wrote in message <hu26a0$hpf$1(a)fred.mathworks.com>... > Mr Stafford, > > Yes, I need f(x) to be the pdf for my generated random variable. Now, maybe I'm confused here, but I don't have that freedom to set p4=p5. Those coefficients come from a nonlinear regression of data, and f(x) fits the data very well. Regarding p1, it is approx equal to = 1/(2*(p3-p2)) * 100. Here are the coeffs: > > p = [4.4782 16.8692 28.0382 6.6695 3.8967]; > > I've plotted f(x) and seem to be always>0. Here's the code > > syms x positive > f = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))); % x is defined from 0 to 80. > norm_c = double(int(f,-1000,1000)); % Normalization coefficient > f = f / norm_c; > F = int(f,x); % integral of f = 1 > ord = subs(F,0); > F = F - ord; > > Given what I wrote above, do I still need the fzero function? I'm stuck at it.. - - - - - - - - - - I don't know what you meant by saying "Regarding p1, it is approx equal to = 1/(2*(p3-p2))*100." When I set p1 equal to 1, and use the values of p2, p3, p4, and p5 you gave, the full integral has a value of 22.338 precisely, and that is also what 2*(p3-p2) is equal to, so you need to set p1 equal to 1/(2*(p3-p2)) = 4.4767e-2 to have a valid density, (as I said before.) Using those values, try plotting f(x) with x = linspace(44,50) and you will see that it goes negative in that range. Even though it does not extend very far into negative values, this creates triple values for the inverse cdf in that region. This would therefore occasionally give wildly inconsistent results from fzero. In answer to your question, the indefinite integral of f(x) has the form: 1/(2*(p3-p2))*(p4*log(cosh((x-p2)/p4))-p5*log(cosh((x-p3)/p5))) and I don't see any way of finding its inverse directly using elementary functions. I think as you say, you are "stuck at it" with fzero or something similar. By the way, I am now convinced 100% that p4 and p5 have to be equal in order to strictly avoid negative values in f(x), as you have defined it. This can be demonstrated rigorously by making use of the expansion of tanh(x) for large values of x: tanh(x) = 1 - 2*exp(-2*x) + 2*exp(-4*x) - 2*exp(-6*x) + ... , or the analogous expansion for large negative x. By going out far enough in the plus or else the minus direction on the x-axis, a negative value for f(x) is certain to be encountered, making a very dubious kind of probability density, at least from a purely mathematical point of view. What it boils down to is that the form you have selected for f(x) cannot be a perfectly valid density with unequal p4 and p5. The true density function would presumably be something a little different, particularly out in the tails. Roger Stafford
From: Roger Stafford on 2 Jun 2010 04:27 I finally got around to working out the details of the cumulative distribution function for your density function f(x). As you recall, your normalized density function was this: f(x) = 1/(2*(p3-p2))*(tanh((x-p2)/p4)-tanh((x-p3)/p5)) where the factor p1 = 1/(2*(p3-p2)) was necessary for the full integral of the density to be one. I have derived three forms for the cdf function which are mathematically identical. F1(x) = 1/2+p6*log(2*cosh((x-p2)/p4))-p7*log(2*cosh((x-p3)/p5)) F2(x) = p6*log(exp(2*(x-p2)/p4)+1)-p7*log(exp(2*(x-p3)/p5)+1) F3(x) = 1+p6*log(1+exp(-2*(x-p2)/p4))-p7*log(1+exp(-2*(x-p3)/p5)) where p6 = p4/(2*(p3-p2)) and p7 = p5/(2*(p3-p2)). In each of them, as x approaches minus infinity, the cdf approaches zero, and as x approaches plus infinity, the cdf approaches one. (As I mentioned earlier, for the values of p2, p3, p4, and p5 you are currently using, F(x) unfortunately goes a small amount beyond one and has to converge on one from above, reflecting a negative density, something self-respecting cdf's are not supposed to do.) The difference between the three formulas is that F2 is designed to handle very large negative values of x, while F3 can handle very large positive values of x. Formula F1(x) could conceivably encounter some loss of accuracy for very large positive or negative values of x as the cosh functions approach infinity, even though the logarithms of these are supposed to cancel out these large values. If you have trouble with F1(x) because of that, you could define an F(x) using F2(x) for, say, negative values of x and F3(x) for positive values. Then you would be assured of not encountering threateningly large intermediate values during the computation. As I have said, I can think of no explicit inverse function for this cdf unless p4 and p5 are equal, so you will have to use something like fzero to solve it for each 'rand' value you generate. If you have trouble with the rate of convergence of fzero, the only suggestion I might make is that you could temporarily use a single intermediate common value in place of p4 and p5, which would then permit an explicit solution for its inverse, with the intent to furnish the required x0 estimate needed in fzero. Then you could hand fzero the version that has p4 and p5 at their original values. I have noticed that if the intermediate value is somewhere between p4 and p5, the two graphs of the respective cdf functions are not terribly far apart, so doing this might speed up the convergence process. Roger Stafford
From: Gonzalo on 2 Jun 2010 07:26
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hu54l2$qhd$1(a)fred.mathworks.com>... > I finally got around to working out the details of the cumulative distribution function for your density function f(x). As you recall, your normalized density function was this: > > f(x) = 1/(2*(p3-p2))*(tanh((x-p2)/p4)-tanh((x-p3)/p5)) > > where the factor p1 = 1/(2*(p3-p2)) was necessary for the full integral of the density to be one. > > I have derived three forms for the cdf function which are mathematically identical. > > F1(x) = 1/2+p6*log(2*cosh((x-p2)/p4))-p7*log(2*cosh((x-p3)/p5)) > F2(x) = p6*log(exp(2*(x-p2)/p4)+1)-p7*log(exp(2*(x-p3)/p5)+1) > F3(x) = 1+p6*log(1+exp(-2*(x-p2)/p4))-p7*log(1+exp(-2*(x-p3)/p5)) > > where p6 = p4/(2*(p3-p2)) and p7 = p5/(2*(p3-p2)). In each of them, as x approaches minus infinity, the cdf approaches zero, and as x approaches plus infinity, the cdf approaches one. (As I mentioned earlier, for the values of p2, p3, p4, and p5 you are currently using, F(x) unfortunately goes a small amount beyond one and has to converge on one from above, reflecting a negative density, something self-respecting cdf's are not supposed to do.) > > The difference between the three formulas is that F2 is designed to handle very large negative values of x, while F3 can handle very large positive values of x. Formula F1(x) could conceivably encounter some loss of accuracy for very large positive or negative values of x as the cosh functions approach infinity, even though the logarithms of these are supposed to cancel out these large values. If you have trouble with F1(x) because of that, you could define an F(x) using F2(x) for, say, negative values of x and F3(x) for positive values. Then you would be assured of not encountering threateningly large intermediate values during the computation. > > As I have said, I can think of no explicit inverse function for this cdf unless p4 and p5 are equal, so you will have to use something like fzero to solve it for each 'rand' value you generate. If you have trouble with the rate of convergence of fzero, the only suggestion I might make is that you could temporarily use a single intermediate common value in place of p4 and p5, which would then permit an explicit solution for its inverse, with the intent to furnish the required x0 estimate needed in fzero. Then you could hand fzero the version that has p4 and p5 at their original values. I have noticed that if the intermediate value is somewhere between p4 and p5, the two graphs of the respective cdf functions are not terribly far apart, so doing this might speed up the convergence process. > > Roger Stafford Mr. Stafford, I'd like to thank you a lot for all the work you did in solving my problem. I couldn't work on it yesterday but I'll go through your notes with care today. Thanks a lot again. I'l prepare a report on the problem and would like to send it to you if you don't mind to your email. Please let me know.. Thanks again, Gonzalo |