Golf

Intro

Minimum and maximum capture velocities

Penner (2002) derives a formula for the maximum capture velocity as: \[v_f < (2R-r)\sqrt{\frac{g}{2R}}\] For a USGA approved golf ball with \(r=2.135cm\) and hole with \(R=5.40cm\) this simplifies to \(v_f<1.31 m/s\).

Rolling resistance

If we assume that the golf ball rolls without slipping or bouncing with initial velocity \(v_i\) conservation of energy gives us: \(\frac{1}{2}mv_i^2 - Fd = \frac{1}{2}mv_f^2\) with constant friction force \(F\) given by \(F = C_rmg\) We can solve for the initial velocity required for \(0 m/s<v_f<1.31 m/s\) as a function of \(d\) Which after some algebra gives: \[\sqrt{2*C_r*9.81}<v_{req}<\sqrt{1.31+2*Cr*9.81*d}\]

x <- 1:25
y1.1 <- sqrt(2*0.05*9.81*x)
y1.2 <- sqrt(2*0.09*9.81*x)
y1.3 <- sqrt(2*0.14*9.81*x)
y2.1 <- sqrt(1.31+2*0.05*9.81*x)
y2.2 <- sqrt(1.31+2*0.09*9.81*x)
y2.3 <- sqrt(1.31+2*0.14*9.81*x)
p <- ggplot() + geom_line(aes(x = x, y = y1.1),  lty = "dashed", col = "green") + 
  geom_line(aes(x = x, y = y1.2), lty = "dashed", col = "blue" ) + 
  geom_line(aes(x = x, y = y1.3), lty = "dashed", col = "red") + 
  geom_line(aes(x = x, y = y2.1), lty = "dashed", col = "green") + 
  geom_line(aes(x = x, y = y2.2), lty = "dashed", col = "blue") +
  geom_line(aes(x = x, y = y2.3), lty = "dashed", col = "red") +
  xlab("x (m)") + ylab("V (m/s)") + ggtitle("Capture Velocity Ranges") + theme_classic()
p

Model 1

data <- read.table("golf_data_new.txt", header=TRUE, skip=2)
data$x <- data$x/3.281 #convert ft to m
data$phat <- data$y/data$n

p <- ggplot(data = data) + geom_point(aes(x = x, y = phat)) + theme_classic()
p

I assume that our golfers are well practiced and will aim to hit the ball at an initial velocity halfway between the two theoretical limits. The velocity they actually hit is modeled as \[v_{target} \sim N(\frac{v_{max}-v_{min}}{2}, \sigma_v)\] Then the probability they will sink the putt (for fixed x) is given as \[Pr(v_{min}< v_{target} < v_{max}) = \Phi(\frac{v_{max}-v_{target}}{\sigma_v})-\Phi(\frac{v_{min}-v_{target}}{\sigma_v})\] where \[v_{min} = \sqrt{2*C_r*9.81*x}\] \[v_{max} = \sqrt{1.31+2*0.09*9.81*d}\] and \[v_{target} = \frac{v_{max} + v_{min}}{2}\] Note that technically, \(d = x - R\) but I don’t think it will make that much difference. Further, I’ll model the \(\sigma_v\) for each bin hierarchically where \[\sigma_{v,i} \sim N(\mu_{sigma_v}, \tau_v)\]. I gave \(\mu_{sigma_v}\) and \(\tau_v\) independent \(N(0,1)\) priors and constrained them to be positive.

J <- length(data$y)
golf.data <- list(x=data$x, y=data$y, n=data$n, J = J)
fit1 <- stan("golf_velocity_1.stan", data = golf.data, chains = 4, warmup = 1000, iter = 5000)

cat(get_stancode(fit1))

## data {
##   int J;
##   int n[J];
##   vector[J] x;
##   int y[J];
## }
## transformed data{
## }
## parameters {
##   real<lower=0> tau_v;
##   real<lower=0> mu_sigma_v;
##   vector<lower=0>[J] sigma_v;
## }
## transformed parameters{
## }
## model {
##   vector[J] p;
##   vector[J] v_max = sqrt(1.31 + 2 * 9.81 * 0.09 * x);
##   vector[J] v_min = sqrt(2 * 9.81 * 0.09 * x);
##   vector[J] v_target = (v_max + v_min)/2;
##   tau_v ~ normal(0,1);
##   mu_sigma_v ~ normal(0,1);
##   sigma_v ~ normal(mu_sigma_v, tau_v);
##   p = Phi((v_max-v_target) ./ sigma_v) - Phi((v_min-v_target) ./ sigma_v);
##   y ~ binomial(n, p);
## }
## generated quantities {
## }

posterior1 <- as.array(fit1)
mcmc_intervals(posterior1[,,1:33])

draws <- extract(fit1)
v_min <- sqrt(2*0.09*9.81*data$x)
v_max <- sqrt(1.31+2*0.09*9.81*data$x)
v_target <- (v_min + v_max)/2
for (i in 1:length(data$x)){
  data$pfit[i] <- mean(pnorm((v_max[i] - v_target[i])/draws$sigma_v[,i]) - pnorm((v_min[i] - v_target[i])/draws$sigma_v[,i]))
}
p <- ggplot(data = data) + geom_point(aes(x = x, y = phat)) + geom_line(aes(x=x, y = pfit)) + theme_classic()
p

sigma_means <- apply(draws$sigma_v, 2, mean)
sigma_95 <- apply(draws$sigma_v, 2, quantile , probs = 0.95)
sigma_05 <- apply(draws$sigma_v,2, quantile, probs = 0.05)
g <- ggplot() + geom_point(aes(x = data$x, y = sigma_means)) + geom_errorbar(aes(x = data$x, ymin = sigma_05, ymax = sigma_95)) +
  geom_abline(slope = 0, intercept = mean(draws$mu_sigma_v), lty = "dashed") + 
  xlab("x(m)") + ylab("Velocity Std Dev") + theme_classic()
g

data$resid <- data$phat-data$pfit
g <- ggplot(data = data) + geom_line(aes(x=x,y=resid)) + geom_abline(slope = 0, intercept = 0, lty = "dotted") + 
  ylab("Residual") + xlab("x (m)") + theme_classic()
g 

Model 2

This time we will estimate \(C_r\) instead of assuming a constant value.

J <- length(data$y)
golf.data <- list(x=data$x, y=data$y, n=data$n, J = J)
fit2 <- stan("golf_velocity_2.stan", data = golf.data, chains = 4, warmup = 1000, iter = 7500, control = list(max_treedepth = 20))

cat(get_stancode(fit2))

## data {
##   int J;
##   int n[J];
##   vector[J] x;
##   int y[J];
## }
## transformed data{
## }
## parameters {
##   real<lower=0> tau_v;
##   real<lower=0> mu_sigma_v;
##   vector<lower=0>[J] sigma_v;
##   real<lower=0> Cr;
## }
## transformed parameters{
## ;
## }
## model {
##   vector[J] p;
##   Cr ~ normal(0,1);
##   tau_v ~ normal(0,1);
##   mu_sigma_v ~ normal(0,1);
##   sigma_v ~ normal(mu_sigma_v, tau_v);
##   vector[J] v_max = sqrt(1.31 + 2 * 9.81 * Cr * x);
##   vector[J] v_min = sqrt(2 * 9.81 * Cr * x);
##   vector[J] v_target = (v_max + v_min)/2;
##   p = Phi((v_max-v_target) ./ sigma_v) - Phi((v_min-v_target) ./ sigma_v);
##   y ~ binomial(n, p);
## }
## generated quantities {
## }

posterior2 <- as.array(fit2)
mcmc_intervals(posterior2[,,1:34])

draws <- extract(fit2)

for (i in 1:length(data$x)){
  v_min <- sqrt(2*draws$Cr*9.81*data$x[i])
  v_max <- sqrt(1.31+2*draws$Cr*9.81*data$x[i])
  v_target <- (v_min + v_max)/2
  data$pfit[i] <- mean(pnorm((v_max - v_target)/draws$sigma_v[,i]) - pnorm((v_min - v_target)/draws$sigma_v[,i]))
}
p <- ggplot(data = data) + geom_point(aes(x = x, y = phat)) + geom_line(aes(x=x, y = pfit)) + theme_classic()
p

data$resid <- data$phat-data$pfit
g <- ggplot(data = data) + geom_line(aes(x=x,y=resid)) + geom_abline(slope = 0, intercept = 0, lty = "dotted") + theme_classic()
g

sigma_means <- apply(draws$sigma_v, 2, mean)
sigma_95 <- apply(draws$sigma_v, 2, quantile , probs = 0.95)
sigma_05 <- apply(draws$sigma_v,2, quantile, probs = 0.05)
g <- ggplot() + geom_point(aes(x = data$x, y = sigma_means)) + geom_errorbar(aes(x = data$x, ymin = sigma_05, ymax = sigma_95)) +
  geom_abline(slope = 0, intercept = mean(draws$mu_sigma_v), lty = "dashed") + 
  xlab("x(m)") + ylab("Velocity Std Dev") + theme_classic()
g

References

Penner, A.J. (2002). “The physics of putting”. Canadian Journal of Physics, 80:1-14.