# Help solving a power equation

Discussion in 'Physics & Math' started by Jennifer Murphy, Apr 8, 2017.

1. ### Jennifer MurphyRegistered Senior Member

Messages:
182
I have a set of N scores, Si, in ascending order. S1 is the smallest, S2 the next smallest, etc.
Code:
   1   2   3   4
S  2   3   5   8

I want to calculate a set of probabilities that are inversely proportional to the scores and that sum to "1". That is, I want S1 to have the highest probability and rest to have probabilities that are proportionately smaller. I accomplished this by first calculating a set of inverse weights, Wi, such that
Code:
Wi = S1 / Si
I then converted those to probabilities, Pi, by
Code:
     1      2      3      4    Sum
S    2      3      5      8
W  1.00   0.67   0.40   0.25   2.32

Then I converted those into probabilities by dividing each one by the sum of the weights, 2.32:
Code:
     1      2      3      4    Sum
S    2      3      5      8
W  1.00   0.67   0.40   0.25   2.32
P  0.43   0.29   0.17   0.11   1.00

Now if P1 (0.43) is greater than some arbitrary value, like 0.30, then I want to set it to that value and scale the other probabilities in proportion to their relative weights, but still have them sum to 1.00.

I can get the result I want by changing the equation for calculating Wi by adding an exponent:
Code:
Wi = (S1 / Si)^F
By trial and error, I found that setting F = 0.26, I get
Code:
     1      2      3      4    Sum
S    2      3      5      8
W  1.00   0.90   0.79   0.70   3.89
P  0.30   0.27   0.23   0.21   1.00

I would like either a way to do that directly or a closed for equation for F. When I try to do that, I get bogged down in complicated (for me) equations involving logs.

I would appreciate any help.

3. ### Jennifer MurphyRegistered Senior Member

Messages:
182
I got this far in the search for a closed form equation for F.

The equation for P1 is:
Code:
P1 = 1 / Sum(Wi)
Sum(Wi) = (S1/S1)^F + (S1/S2)^F + (S1/S3)^F
P1 = 1 / ( (S1/S1)^F + (S1/S2)^F + (S1/S3)^F )
(S1/S1)^F + (S1/S2)^F + (S1/S3)^F  = 1/P1

I have come up short for a way to solve this for F. Normally, when trying to solve for an exponent, I would take the logs of both sides, but I don;t know how to get F out of that log of sums.

Can anyone suggest a way to solve this equation for F ?

5. ### rpennerFully WiredValued Senior Member

Messages:
4,833
So you have a vector of positive numbers S and you want a second vector of positive numbers P and an positive exponent K with conditions:
max(P_i) = P0
sum(P_i) = 1
S_i P_i^K = S_j P_j^K

This will not always be possible. If S is constant, P0 has to be a particular value.

If P_i = W S_i^-K then
S_i P_i^K = W = S_j P_j^K
max(P_i) = W min(S_i)^-K

Thus W(K) = P0 min(S_i)^K
Thus sum(P_i) = P0 min(S_i)^K sum(S_i^-K)

So we want to solve for K:
1 = P0 min(S_i)^K sum(S_i^-K)​
for P0 = 3/10

This in general is a transcendental equation so you will want a general solver and contingency for when there is no solution or no positive solution. I get K ≈ 0.285661799704292683848547628090811554694676399989531... ≈ 2/7

Using K=2/7 (which is just and approximation), you get W = (3/10) 2^(2/7) and (3/10) 2^(2/7) ( 2^(-2/7) + 3^(-2/7) + 5^(-2/7) + 8^(-2/7) ) = 0.9999685...

So using Mathematica, we can get a highly precise solution:

Code:
S = { 2, 3, 5, 8 };
P0 = 3/10 ;
P = P0 * Min[S]^x * ( # ^ -x ) & /@ S ;
ans = Solve[ Total[ P ] == 1, x, Reals] ;
N[P , 79] /. ans[[1]]

Code:
{0.3000000000000000000000000000000000000000000000000000000000000000000000000000000,
0.267189082563102067083882226912844813277540705949340595586058500878920594709311,
0.230911198638986395604416849051029728933476685369502782406193969753574064993534,
0.201899718797911537311700924036125457788982608681156622007747529367505340297155}
Replacing some of the scary syntax with the long forms gives:
Code:
S = { 2, 3, 5, 8 };
P0 = 3/10 ;
P = P0 * Min[S]^x * Map[ Function[ s, s ^ -x ] , S] ;
ans = Solve[ Total[ P ] == 1, x, Reals] ;
ReplaceAll[N[P , 79] , Part[ans,1] ]


Last edited: Apr 13, 2017