data=read.csv("karpur.csv")
head(data)
## depth caliper ind.deep ind.med gamma phi.N R.deep R.med SP
## 1 5667.0 8.685 618.005 569.781 98.823 0.410 1.618 1.755 -56.587
## 2 5667.5 8.686 497.547 419.494 90.640 0.307 2.010 2.384 -61.916
## 3 5668.0 8.686 384.935 300.155 78.087 0.203 2.598 3.332 -55.861
## 4 5668.5 8.686 278.324 205.224 66.232 0.119 3.593 4.873 -41.860
## 5 5669.0 8.686 183.743 131.155 59.807 0.069 5.442 7.625 -34.934
## 6 5669.5 8.686 109.512 75.633 57.109 0.048 9.131 13.222 -39.769
## density.corr density phi.core k.core Facies phi.core.frac
## 1 -0.033 2.205 33.9000 2442.590 F1 0.339000
## 2 -0.067 2.040 33.4131 3006.989 F1 0.334131
## 3 -0.064 1.888 33.1000 3370.000 F1 0.331000
## 4 -0.053 1.794 34.9000 2270.000 F1 0.349000
## 5 -0.054 1.758 35.0644 2530.758 F1 0.350644
## 6 -0.058 1.759 35.3152 2928.314 F1 0.353152
data = data[order(data$k.core, decreasing=TRUE), ]
k = data$k.core
sample = c(1: length(k))
k_percent = (sample * 100 / length(k))
xlab = "Portion of Total Samples Having Larger or Equal K"
ylab = "Permeability k (md)"
plot(k_percent, k, log = 'y', xlab = xlab, ylab = ylab, pch = 10, cex = 0.5, col = "#16423C")

log.k = log(k)
model = lm(log.k ~ k_percent)
plot(k_percent, log.k, xlab = xlab , ylab = ylab, pch = 10, cex = 0.5, col = "#16423C")
abline(model, col = "red", lwd = 2)

N_data = data.frame(k_percent = c(50, 84.1))
predict_values = predict(model, N_data)
HI = (predict_values[1] - predict_values[2]) / predict_values[1]
HI
## 1
## 0.2035464