From 42aed2eec12d2dc6f3cd976f1851cb3094053924 Mon Sep 17 00:00:00 2001
From: field <field@igskci164lm027.myverizon>
Date: Wed, 14 Aug 2024 13:25:37 -0700
Subject: [PATCH] Added a model and fixed a bug in the Lognormal distribution
 creation

---
 .../kevin/ucerf3/eal/LossCOV_Model.java       | 35 +++++++++++++++++--
 1 file changed, 32 insertions(+), 3 deletions(-)

diff --git a/src/main/java/scratch/kevin/ucerf3/eal/LossCOV_Model.java b/src/main/java/scratch/kevin/ucerf3/eal/LossCOV_Model.java
index 9cc3a555..3e2091b7 100644
--- a/src/main/java/scratch/kevin/ucerf3/eal/LossCOV_Model.java
+++ b/src/main/java/scratch/kevin/ucerf3/eal/LossCOV_Model.java
@@ -3,6 +3,7 @@
 import java.awt.Color;
 
 import org.apache.commons.math3.distribution.LogNormalDistribution;
+import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
 import org.jfree.data.Range;
 import org.opensha.commons.data.function.ArbitrarilyDiscretizedFunc;
 import org.opensha.commons.data.function.DiscretizedFunc;
@@ -15,18 +16,36 @@
 
 public enum LossCOV_Model {
 	
-	PORTER_POWER_LAW_2020_09_01 {
+//	// not implementing this because it requires the x-axis values of the func supplied in calcLossExceedanceProbs to be modified too
+//	PORTER_POWER_LAW_2024_08_09 { // meanLoss is loss (in dollars) divided by total portfolio value in units of $1000
+//		@Override
+//		public double getCOV(double meanLoss) {
+//			return Math.min(2.026, 0.9832*Math.pow(meanLoss, -0.117));
+//		}
+//	},
+	PORTER_POWER_LAW_2020_09_01 { // Here, meanLoss is loss in units of $1000
 		@Override
 		public double getCOV(double meanLoss) {
 			return 4.4751*Math.pow(Math.max(1000, meanLoss), -0.1152);
 		}
+	},
+	PORTER_POWER_LAW_2020_09_01_fixed { // fixed to be consistent with Keith's manuscripts; Here, meanLoss is loss in units of $1000
+		@Override
+		public double getCOV(double meanLoss) {
+			return 4.546*Math.pow(Math.max(1000, meanLoss), -0.117);
+		}
 	};
+
 	
 	public abstract double getCOV(double meanLoss);
 	
 	public LogNormalDistribution getDistribution(double meanLoss) {
 		Preconditions.checkState(meanLoss > 0d);
-		return new LogNormalDistribution(Math.log(meanLoss), getCOV(meanLoss));
+		double cov = getCOV(meanLoss);
+		double sigma = Math.sqrt(Math.log(cov*cov+1));
+		double mu = Math.log(meanLoss)-(sigma*sigma/2);
+
+		return new LogNormalDistribution(mu, sigma);
 	}
 	
 	public DiscretizedFunc calcLossExceedanceProbs(DiscretizedFunc xVals, double meanLoss) {
@@ -46,7 +65,17 @@ public double calcLossExceedanceProb(double x, double meanLoss) {
 	}
 	
 	public static void main(String[] args) {
-		LossCOV_Model model = PORTER_POWER_LAW_2020_09_01;
+		LossCOV_Model model = PORTER_POWER_LAW_2020_09_01_fixed;
+		
+		// test lognormal distribution
+        double mean = 1e5;
+        double cov = model.getCOV(mean);
+        double samples[] =      model.getDistribution(mean).sample(100000000);
+        DescriptiveStatistics stats = new DescriptiveStatistics(samples);
+        double meanFromDist = stats.getMean();
+        double covFromDist = stats.getStandardDeviation()/meanFromDist;
+       System.out.println("mean="+mean+";\tmeanFromDist="+meanFromDist+";\tcov="+cov+";\tcovFromDist="+
+    		   covFromDist+"; meanRatio="+(float)(meanFromDist/mean)+"; covRatio="+(float)(covFromDist/cov));
 		
 		EvenlyDiscretizedFunc logFunc = new EvenlyDiscretizedFunc(0d, 8d, 500);
 		ArbitrarilyDiscretizedFunc func = new ArbitrarilyDiscretizedFunc();