Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Generate random numbers from a particular function

Subject: Generate random numbers from a particular function

From: Gonzalo

Date: 31 May, 2010 22:10:19

Message: 1 of 10

Hi all,

I need to generate RN's from a double hyperbolic tangent function.

f(x) = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))

I can't do it using the inverse-transform method because it's not possible to solve for x. Does anybody know of a routine that works the Composition, or convolution or any other methods??

Thanks

Subject: Generate random numbers from a particular function

From: Walter Roberson

Date: 31 May, 2010 23:33:17

Message: 2 of 10

Gonzalo wrote:

> I need to generate RN's from a double hyperbolic tangent function.
>
> f(x) = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))
>
> I can't do it using the inverse-transform method because it's not
> possible to solve for x. Does anybody know of a routine that works the
> Composition, or convolution or any other methods??

x can be solved for numerically given the other parameters. The key value to
be solved for is,

RootOf(2*_Z*p(5)-2*p(3)+2*p(2)-p(4)*ln((p(1)*exp(_Z)^2+p(1)+exp(_Z)^2-1)/(p(1)*
exp(_Z)^2+p(1)-exp(_Z)^2+1)))

where RootOf is a notation indicating that the value _Z should be found such
that the expression evaluates to 0 at _Z .

However, for some combinations of parameters, some of the x might be
imaginary. For example, p(1)=1/2, p(2)=1/3, p(3)=1/5, p(4)=1/7 and p(5) from
about 0.18 to about 0.54, whereas with p(1)=2, p(2)=3, p(3)=5, p(4)=7 then in
my experiments I do not see any p(5) that would make the expression imaginary.

Subject: Generate random numbers from a particular function

From: Gonzalo

Date: 1 Jun, 2010 00:19:05

Message: 3 of 10

Thanks Walter,

I found the coefficients with a nonlinear regression to be:

p = [4.4786 16.8696 28.0379 3.3347 1.9488];

Could you explain a little more how to solve numerically??

thanks a lot,

Gonzalo

Walter Roberson <roberson@hushmail.com> wrote in message <hu1h27$e6l$1@canopus.cc.umanitoba.ca>...
> Gonzalo wrote:
>
> > I need to generate RN's from a double hyperbolic tangent function.
> >
> > f(x) = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))
> >
> > I can't do it using the inverse-transform method because it's not
> > possible to solve for x. Does anybody know of a routine that works the
> > Composition, or convolution or any other methods??
>
> x can be solved for numerically given the other parameters. The key value to
> be solved for is,
>
> RootOf(2*_Z*p(5)-2*p(3)+2*p(2)-p(4)*ln((p(1)*exp(_Z)^2+p(1)+exp(_Z)^2-1)/(p(1)*
> exp(_Z)^2+p(1)-exp(_Z)^2+1)))
>
> where RootOf is a notation indicating that the value _Z should be found such
> that the expression evaluates to 0 at _Z .
>
> However, for some combinations of parameters, some of the x might be
> imaginary. For example, p(1)=1/2, p(2)=1/3, p(3)=1/5, p(4)=1/7 and p(5) from
> about 0.18 to about 0.54, whereas with p(1)=2, p(2)=3, p(3)=5, p(4)=7 then in
> my experiments I do not see any p(5) that would make the expression imaginary.

Subject: Generate random numbers from a particular function

From: Roger Stafford

Date: 1 Jun, 2010 00:44:04

Message: 4 of 10

"Gonzalo " <glpita@gmail.com> wrote in message <hu1c4b$5t7$1@fred.mathworks.com>...
> Hi all,
>
> I need to generate RN's from a double hyperbolic tangent function.
>
> f(x) = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))
>
> I can't do it using the inverse-transform method because it's not possible to solve for x. Does anybody know of a routine that works the Composition, or convolution or any other methods??
>
> Thanks

  I am guessing that you mean for f(x) to be the probability density function for your random variable, though you didn't say so. Is that correct? If so, the integral of f(x) will be a linear combination of two log cosh terms, and for that reason the inverse of the cdf should be comparatively easy to obtain from fzero. If this fzero process isn't too slow for you, you could then generate the random variable using the rand function in connection with fzero in accordance with the standard method.

  Of course your five p-values are subject to some constraints. The integral of f(x) must be unity, and there should be no x for which f(x) becomes negative. You would have to be sure of these conditions before attempting to generate a valid random variable.

Roger Stafford

Subject: Generate random numbers from a particular function

From: Gonzalo

Date: 1 Jun, 2010 01:33:04

Message: 5 of 10

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hu1l4k$ob7$1@fred.mathworks.com>...
> "Gonzalo " <glpita@gmail.com> wrote in message <hu1c4b$5t7$1@fred.mathworks.com>...
> > Hi all,
> >
> > I need to generate RN's from a double hyperbolic tangent function.
> >
> > f(x) = p(1) .* (tanh((x - p(2)) ./ p(4))- tanh((x - p(3)) ./ p(5))
> >
> > I can't do it using the inverse-transform method because it's not possible to solve for x. Does anybody know of a routine that works the Composition, or convolution or any other methods??
> >
> > Thanks
>
> I am guessing that you mean for f(x) to be the probability density function for your random variable, though you didn't say so. Is that correct? If so, the integral of f(x) will be a linear combination of two log cosh terms, and for that reason the inverse of the cdf should be comparatively easy to obtain from fzero. If this fzero process isn't too slow for you, you could then generate the random variable using the rand function in connection with fzero in accordance with the standard method.
>
> Of course your five p-values are subject to some constraints. The integral of f(x) must be unity, and there should be no x for which f(x) becomes negative. You would have to be sure of these conditions before attempting to generate a valid random variable.
>
> Roger Stafford

Roger,

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)??

thanks

Subject: Generate random numbers from a particular function

From: Roger Stafford

Date: 1 Jun, 2010 04:36:03

Message: 6 of 10

"Gonzalo " <glpita@gmail.com> wrote in message <hu1o0g$44h$1@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

Subject: Generate random numbers from a particular function

From: Gonzalo

Date: 1 Jun, 2010 05:37:04

Message: 7 of 10

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hu22nj$jv4$1@fred.mathworks.com>...
> "Gonzalo " <glpita@gmail.com> wrote in message <hu1o0g$44h$1@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..

Subject: Generate random numbers from a particular function

From: Roger Stafford

Date: 1 Jun, 2010 09:27:03

Message: 8 of 10

"Gonzalo " <glpita@gmail.com> wrote in message <hu26a0$hpf$1@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

Subject: Generate random numbers from a particular function

From: Roger Stafford

Date: 2 Jun, 2010 08:27:14

Message: 9 of 10

  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

Subject: Generate random numbers from a particular function

From: Gonzalo

Date: 2 Jun, 2010 11:26:04

Message: 10 of 10

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hu54l2$qhd$1@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

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us