Help solving a power equation

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

  1. Jennifer Murphy Registered Senior Member

    Messages:
    175
    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.
     
  2. Google AdSense Guest Advertisement



    to hide all adverts.
  3. Jennifer Murphy Registered Senior Member

    Messages:
    175
    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 ?
     
  4. Google AdSense Guest Advertisement



    to hide all adverts.
  5. rpenner Fully Wired Staff 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
  6. Google AdSense Guest Advertisement



    to hide all adverts.

Share This Page