Means of Getting Means

Finding the average of a set of values (of samples) is something we all know. Annie has 50 dollars, Becky has 30, and Carol has 10. What is the average amount of money they hold? You sum the numbers and divide by the number of sums, i.e., (50+30+10)/3 = 90/3 = 30. This is mean, but more specifically, the arithmetic mean.

\bar x = \frac{\sum x}{n}

Now look at this: I have invested 100 dollars in the market. It gave 10% in the first year, 20% in the second, and 30% in the third. What is the average annualised rate of return? The first instinct is to take the arithmetic mean of the three numbers and report it as 20%. Well, is that correct? Let’s check.

Money at the end of year 1 = 100 x (1.1) = 110
Money at the end of year 2 = 110 x (1.2) = 132
Money at the end of year 3 = 132 x (1.3) = 171.6

Now, let’s calculate the final amount using the proposed mean, i.e., 20%. The calculation led to 100x(1.2)3 = 172.8; not quite the same!

Here, we estimate the geometric mean by multiplying the returns and taking the cube root. Mean = [(1.1)x(1.2)x(1.3)]1/3 = 1.7163 = 1.197216; the average return is 0.1972 or 19.72%.
The average return using the new (geometric) mean is 100x(1.197216)3 = 171.6.

\sqrt[n]{\Pi x}

Abby ran the first lap at 10 km/hr speed and the second at 12 km/hr. What is Abby’s average speed? The answer is neither (10+12)/2 nor (10×12)1/2. It is 2/(1/10 + 1/12) = 10.9. It is known as the harmonic mean and is given by:

\frac{n}{\sum(1/x)}

Means of Getting Means Read More »

Coefficient of Variation

We know standard deviation quantifies the spread of data from the mean. It is an absolute measure, and therefore, it has the same unit as the variable. However, it is limited when comparing two standard deviations of samples having different means. For example, compare the spread of the following two samples: the spending habits of two types of households.

Sample 1 (high income = mean income 400,000): standard deviation = 40,000
Sample 2 (low income = mean income 10,000): standard deviation = 1,000

It appears that the variability of data in the high-income category is higher. Or is it?

In such cases, non-dimensionalizing or standardization of the spread comes in handy. The Coefficient of Variation (CV) is one such. It is the ratio between the standard deviation and the mean. Since both possess the same unit, the ratio is dimensionless.

CV = Standard Deviation / Mean

Let’s calculate the CVs of sample 1 and sample 2.
CV1 (high income) = 400,000 / 40,000 = 10
CV2 (low income) = 10,000 / 1,000 = 10

Practically the same relative variability.

Coefficient of Variation Read More »

Studies Everyday

Anna must select seven classes out of the 30 non-overlapping courses. Each course happens once a week and is distributed equally from Monday through Friday (6 courses a week). What is the probability that she has courses every day (Monday through Friday)?

\frac{(_5C_2)(_6C_2)^2 6^3 + (_5C_1)(_6C_3) 6^4}{_{30}C_7} = \frac{114}{377} = 0.302

Studies Everyday Read More »

Probability of Circle in a Square

The probability of hitting inside the circle inscribed in a square is a well-known problem. Looking at the geometries, one can find that the required probability is pi/4.

If X is the radius of the circle, then the diameter, 2X, is equivalent to the side of the square. The required probability is the area of the circle (green) over the area of the square. I.e., pi x X2 / (2X) x (2X) = pi / 4. Or four times this probability is pi.

We will do an R code to simulate this case.
Step 1: Create random uniformly distributed points between two equal intervals (-1, 1) to produce a 2 x 2 square.

pts <- 10000
xx <- runif(pts, -1, 1)
yy <- runif(pts, -1, 1)
plot(xx, yy, asp = 1, col = "blue")

Step 2: Select only those points that satisfy the following conditions for a circle: X2 + Y2 < 1, 1 being the circle’s radius that fits inside the square.

for(i in 1:pts){
   if(xx[i]^2 < 1 - yy[i]^2) {
      xx1[i] <- xx[i]
   }else{
      xx1[i] <- 0    
   }
 }
yy1[which(xx1 == 0)] <- 0
plot(xx1, yy1, asp = 1, type = "p", col = "red")

Combining two plots,

The ratio of the number of points inside the circle to that of the square multiplied by 4 is:

4*length(xx1[abs(xx1) > 0 ]) / length(xx)
3.1456

Tailpiece

Increase the number of points to 100,000, and you get a much cleaner picture:

Probability of Circle in a Square Read More »

Romeo, Juliet and Cofee Meet

Romeo and Juliet plan to meet over coffee between 3 PM and 4 PM. Not known for punctuality, either is willing to wait 15 minutes for the other to show up. What is the probability they will have coffee together?

We will do a graphical solution followed by simulations.

The required probability is the chance of two people meeting inside the shaded area of the circle. It is estimated as 1 – chance of both not making it inside the region = 1 – (3/4)x(3/4) = (16-9)/16 = 7/16 = 0.437.

R Simulations

itr <- 1000000
time_band <- 15

  romeo <- runif(itr, 0, 60)
  juliet <- runif(itr, 0, 60)  
  mean(abs(romeo - juliet) <= time_band)
0.437

Romeo, Juliet and Cofee Meet Read More »

One Crime and Two Suspects

Police know that the crime was committed by one of the two suspects, A and B, with equal evidence against them initially. Further investigation suggests that the blood type of the guilty is found in 1/10 of the population. If A’s blood matches this type, whereas B’s is unknown, what is the probability that A is guilty?

The probability that A is guilty, given his blood matches = P(A|M)
P(A|M) = P(M|A) x P(A) /(P(M|A) x P(A) + P(M|B) x P(B))
P(A|M) = 1 x (1/2) /[1 x (1/2) + (1/10)x (1/2)] = 10/11

One Crime and Two Suspects Read More »

Story Proof for Vandermonde’s identity

Vandermonde’s identity is a famous relationship between binomial coefficients. It has the following form:

_{m+n}C_k = \sum\limits_{j=0}^k _{m}C_j * _{n}C_{k-j}

Blitzstein and Hwang, in their book Introduction to Probability, give a story proof to the identity. Suppose m boys and n girls want to form a committee of k members; there are (m+n)Ck possibilities for doing it. It forms the left-hand side of the identity.

For every j number of boys in the committee, there will be (k-j) girls in it. The number of possibilities is mCj x nCk-j. The required possibility will be the sum of all those possibilities where j ranges from 0 (no boys) to k (only boys).

Story Proof for Vandermonde’s identity Read More »

Bayesian Medicine

A new medicine is developed against a disease to improve the incumbent treatment that is effective in 50% of patients. The expectation is that there is a 2/3 chance that the new treatment will be effective on 60% of patients and a 1/3 chance that the new treatment will be as good as the older one (effective on 50% of patients). The new treatment was given to 20 random patients in the first field trial and was found effective for 15.

With the above information, what is the probability that the new treatment is better than the existing one?

The objective is to find the probability that the new medicine (N) is better, given it was effective on 15 patients in the trial. The notation for this is P(N|X=15) and using Bayes’ theorem,

P(N|X=15) = P(X=15|N) x P(N) / [P(X=15|N) x P(N) + P(X=15|Nc) x P(Nc)]

P(X=15|N) = Probability of getting 15 successes given the probability is 0.6 = 20C15 x (0.6)15 x (0.4)5 (Binomial theorem)
P(X=15|Nc) = Probability of getting 15 successes given the probability is 0.5 = 20C15 x (0.5)15 x (0.5)5 (Binomial theorem)
P(N) = 2/3
P(Nc) = 1/3

dbinom(15, 20, prob = 0.6)*(2/3) / (dbinom(15, 20, prob = 0.6)*(2/3) + dbinom(15, 20, prob = 0.5)*(1/3))
0.91

Bayesian Medicine Read More »

Linear Programming in R – Integer Value Optimisation

In the earlier linear optimization problem, we got the solution as X1 (product 1) = 24 and X2 (product 2) = 16. The interesting thing to notice here is that both are whole numbers; you can’t make and sell fractional products! It may only be the case in some cases. Here is an example.

A factory makes three products: chairs, tables and sofas.
One chair makes $150 profit, requires 7 hours to produce, and 3 m2 of space to store
One table makes $250 profit, requires 9 hours to produce, and 6 m2 of space to store
One sofa makes $300 profit, requires 11 hours to produce, and 7 m2 of space to store

Now, the constraints. The resources available for the factory are only sufficient to cover up to 80 hours a week. The warehouse has an area of 40 m3, and based on historical data, a maximum of 4 chairs sell in a week. Estimate an optimised production plan.

Here is an R code and the results.

f.obj <- c(150, 250, 300)
f.con <- matrix (c(7, 9, 11, 3, 6, 7, 4, 0, 0), nrow = 3, byrow = TRUE)
f.dir <- c("<=", "<=", "<=")
f.rhs <- c(80, 40, 8)
results <- lp ("max", f.obj, f.con, f.dir, f.rhs)
results
results$solution
Success: the objective function is 1757.143 
[1] 2.000000 0.000000 4.857143

While making two chairs a week seems reasonable, it doesn’t make any sense to make 4.85 sofas and expect someone to buy a 0.85 item. At the same time, rounding off to 5 will end up violating the two constraints (on resource and space).
So, we will specify integer constraints to variables 1, 2 and 3 by adding ‘int.vec = ‘ in the options.

f.obj <- c(150, 250, 300)
f.con <- matrix (c(7, 9, 11, 3, 6, 7, 4, 0, 0), nrow = 3, byrow = TRUE)
f.dir <- c("<=", "<=", "<=")
f.rhs <- c(80, 40, 8)
results <- lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec = c(1,2,3))
results
results$solution
Success: the objective function is 1750 
[1] 2 1 4

Reference

Linear Optimization: Desmond C. Ong

Linear Programming in R – Integer Value Optimisation Read More »

Linear Programming in R – Shadow Prices

In the last post, we have seen how R performs linear programming and estimates the limits of the coefficients of the objective function over which the optimal solution remains the same. That was one kind of sensitivity check. The second way is to vary the constraints and observe how the objective function varies. It is the shadow price (of a constraint) and is defined as the change in the objective function value for a unit increase on the right-hand side of the constraint.

To recap, the following is the summary of the objective function and the constraints.

Maximise 45 X1 + 50 X2 (Objective function)
500 X1 + 500 X2 </= 20000 (constraint 1)
750 X1 + 625 X2 </= 42000 (constraint 2)
150 X1 + 100 X2 </= 10400 (constraint 3)
200 X1 + 300 X2 </= 9600 (constraint 4)
X1  </= 50 (constraint 5)
X2  </= 20 (constraint 6)
X1  >/= 0 (constraint 7)
X2  >/= 0 (constraint 8)

Also, recall that the above set of equations gave the maximum profit (the objective function value) of 1880, and the optimal production quantities are alpha = 24 and beta = 16.

The shadow price of constraint 1 is the change in maximum profit from 1880 when the right-hand side of constraint 1 is increased from 20000 to 20001. We can run the code and find out what it is.

f.obj <- c(45, 50)
f.con <- matrix (c(500, 500, 750, 625, 150, 100, 200, 300, 1, 0, 0, 1), nrow = 6, byrow = TRUE)
f.dir <- c("<=", "<=")
f.rhs <- c(20001, 42000, 10400, 9600, 50, 30)

results <- lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)
results
Success: the objective function is 1880.07

The value increased by 0.07 (1880.07 – 1880.0). The R can give all the shadow values using the following command. Note that you have added the attribute, ‘compute.sens=TRUE’, in the earlier code.

results$duals
 0.07 0.00 0.00 0.05 0.00 0.00 0.00 0.00

The numbers are the shadow values from constraint 1 to constraint 8. The results show that only a change of one unit of constraint 1 and constraint 4 makes changes in the objective function value. They are, respectively, supply of ingredient A and supply of ingredient D. The former has a shadow value of 0.07, and the latter has 0.05.

Linear Programming in R – Shadow Prices Read More »