Solving a non-linear system of equations. (1 Viewer)

[Lazarus]

Can one of you help me out with this please? turtle? drbuchanan?


I have a set of figures (more than 10) which all fall within the range 0 <= w <= 1. I have mapped the figures into the range (-infinity, infinity) by the transform:

X = ln[(w + 0.0001)/(1 - w + 0.0001)]

I wish to use the inverse transformation with values a and b chosen such that the values returned by

Z = 50(e^{aX + b})/(1 + e^{aX + b})

have a specified mean (M) and standard deviation (S).

It is possible to form two equations:

1) M = 1/n * (Z₁ + Z₂ + Z₃ + ... Z_n)

2) S² = 1/n * [(Z₁ - M)² + (Z₂ - M)² + (Z₃ - M)² + ... + (Z_n - M)²]

Substitute the expression for Z given earlier, and we have a non-linear system of equations.

How does one solve such a system for a and b? I am told I can use numerical iterative procedures (e.g. Newton-Raphson), but I'm not really sure how to do this, and most of the websites I've found on the topic are either too simple or too complex.

Any help would be greatly appreciated.

 

[KeypadSDM]

Originally posted by Lazarus
I am told I can use numerical iterative procedures (e.g. Newton-Raphson)

The only thing in there I actually understood. You can trisect the angle by iterative methods

Why are you posting in the 4 unit board tho?

 

[Lazarus]

Ergh, did I not explain it very well? Sorry, heh. Which parts are difficult to follow?

Originally posted by KeypadSDM
Why are you posting in the 4 unit board tho?

D'you think it's more likely to receive an answer in another forum?

 

[turtle_2468]

Originally posted by Lazarus
Can one of you help me out with this please?

I'll try But first a question:

by the transform:

X = ln[(w + 0.0001)/(1 - w + 0.0001)]
do you mean
X=ln[(w+0.0001)/(1-(w+0.0001)]?

 

[turtle_2468]

Okay. It is solvable using iterative methods for 1 variable... (when I say lots of rearranging, it works. If you want the proof, say if this result is useful, ask me and I will put it up. Just lots of algebra..)
after lots of rearranging, letting m_i denote 1/(1+e^(ax_i+b)) for convenience, given n, s and m you can get
sum of m_i = k_1
sum of (m_i)^2=k_2
where k_1, k_2 are constants.
Then suppose you fix a. Let e^b=c. Then the equation becomes one purely in variable c for say the first equation, and you can solve for c using iterative methods.
NB When it says newton-raphson, it just means the plain vanilla 4U newton's method.

PS also note that if you get the z substitution in the form or something like 50(1+e^(ax+b))/(e^(ax+b)) instead of the other way around it might make the algebra a bit easier... because right now I'm stuck with 1/(1+ce^(x_i)), and if that happened (I don't know whether you can fiddle with the first substitution getting X so you can do that) then the bottom would be without the 1, always nice!

I'll try to see whether any nice substitutions give a good way to solve for both b and c (hopes are slim!) and get back hopefully before tonight.

 

[Lazarus]

Thanks turtle.

Originally posted by turtle_2468
by the transform:

X = ln[(w + 0.0001)/(1 - w + 0.0001)]

do you mean

X=ln[(w+0.0001)/(1-(w+0.0001)] ?

Aren't they the same? The numerator of the fraction within the logarithm should be (w + 0.0001) and the denominator should be (1 - w + 0.0001).

Then suppose you fix a. Let e^b=c. Then the equation becomes one purely in variable c for say the first equation, and you can solve for c using iterative methods.

Okay, I followed everything up to here... I can see the rearranging you've done, etc. The problem is with fixing a - what do we fix it at, and how can we know it will provide a solution (regardless of the value of b)? Or perhaps it might work for only one of the equations; i.e. fixing a will allow a solution to be found for the M equation, but the solution turns out to be different for the S equation (and they must be the same). Hrmmm. Well, you've made progress at least, and given me some more to chew on. Thanks.

NB When it says newton-raphson, it just means the plain vanilla 4U newton's method.

I became confused when trying to extend it to equations in two unknowns... it involved solving a Jacobian matrix (!). Not really sure how they differ to normal matrices, but they sound harder.

 

[turtle_2468]

I meant ln[(w+0.0001)/(1-w-0.0001)].
edit: I guess it doesn't make the thing any nicer after all... although you don't get stuck with an unnice 1.0002 in the numerator at the end
eg
(w+0.0001)/(1-w-0.0001)=(1/(0.999-w))-1
But then that can just be factored out in the final thing in the other case so doesn't matter as much...

 

[Lazarus]

Oh, right, I see - no.

 

[turtle_2468]

Eww. Jacobian matrices are not nice, you encounter them in several variable. I've done most of that course (and I definitely remember doing the bit about jacobians), can you send me the link with newton & relation to jacobian? Then I can dust off my books and get back to work
well, not that I have any books, but notes scrounged from maths will suffice

 

[Lazarus]

Check out this small excerpt I found:
http://www.amath.washington.edu/courses/301-autumn-2003/301lec17.pdf

Specifically, from the end of the fourth page through to the fifth.

 

[turtle_2468]

think I understood the pdf file, but packing to go to blue mountains tomorrow for a week, will relook up jacobians when I get back... (26th) give me a virtual kick if I haven't replied by say 28th or so.

 

[maniacguy]

Quick comments, having just seen the problem 10 minutes ago (sorry, I haven't visited this forum much in recent times. My bad!):

(1) David, he means +0.0001 in both numerator and denominator (this is just to make sure that you can't get 0 or a negative quantity in there for 0 <= w <= 1)

(2) Given the reason why you would have that small 0.0001 it may actually be easier to treat the problem simply by ignoring it, since it shouldn't make a huge difference so that:
0 <= w <= 1
X = ln(w/(1-w))
Y = exp(a*X + b)
Z = 50*Y/(1+Y)

E(Z) = M
Var(Z) = S^2

Thus, treating it as a statistics question, we can easily get it to the point (well, without too much trouble, anyhow) where:

T = A*[w/(1-w)]^B
R = 1/(1+T)

E(R), Var(R) are known constants

Now solve for A, B and you're home (if the w/(1-w) makes you really uncomfortable, then congratulations, you're on your way to becoming a pure mathematician. Turn it into something like (w+delta)/(1-w+delta) as delta -> 0 if you want)

I'm fairly sure there are statistical methods for approximating the constants in this case, provided we can postulate w as representing a known distribution (in this case, the empirical distribution function produced by the sample values). If not, then I'm sure some of the applied mathematicians can help out!

Hope that helps go in vaguely a relevant direction (e.g. if you have friends who have been doing the approximation kind of thing for a while, or are proficient in a computer program. You could probably write a C program to simulate for a while!!), and I'll keep thinking a bit more about it.

edit:
Thinking of this as a random variable situation - and transforming that variable - could also help with the Jacobians, used in such transformations...

 

[KeypadSDM]

Man, I'll be thinking like this in a year. Hopefully.

 

[turtle_2468]

try 2 years if you're talking about maniacguy..

 

[Dash]

Man, thats hectic stuff...
I can't really understand much

 

[maniacguy]

A further comment on what I said about approximation:

If you give w a known density (based on what you've said about all the figures falling between 0 and 1, without much else, I'm going to choose the Uniform(0,1) density), then you can get:

INTEGRAL[w=0..1, 1/(1+A*(w/(1-w))^B) dw] is known
INTEGRAL[w=0..1, 1/(1+A*(w/(1-w))^B)^2 dw] is known

And most decent mathematical software will let you try and approximate this (or even evaluate it as a function of A and B if it's really good!)

you can also simplify the above by setting m = w/(1-w) and finding the density of m (this is straightforward, so I'm leaving it as an exercise aka I'm too lazy to do it myself right now)

 

[Lazarus]

Thanks for your help, maniacguy.

I managed to get some good approximations with Matlab - they seem to work!

 

[Grey Council]

o shit. Do we have to be as good/as naturally talented at mathematics as Maniacguy to do actuarial studies? day-yum, i was thinking of doing actuarial studies in uni.

i'm scared.
*mweep*

 

[turtle_2468]

Umm... you don't

 

[Grey Council]

good. so i have a hope after all