时间单位线计算优化
cxw
parent
61e32a07bf
commit
6d71b6d03b
|
|
@ -1,12 +1,16 @@
|
|||
package com.gunshi.project.xyt.controller;
|
||||
|
||||
import com.baomidou.mybatisplus.core.conditions.query.QueryWrapper;
|
||||
import com.baomidou.mybatisplus.core.conditions.update.UpdateWrapper;
|
||||
import com.baomidou.mybatisplus.core.toolkit.CollectionUtils;
|
||||
import com.baomidou.mybatisplus.core.toolkit.IdWorker;
|
||||
import com.baomidou.mybatisplus.core.toolkit.ObjectUtils;
|
||||
import com.baomidou.mybatisplus.core.toolkit.StringUtils;
|
||||
import com.baomidou.mybatisplus.extension.plugins.pagination.Page;
|
||||
import com.gunshi.core.result.R;
|
||||
import com.gunshi.project.xyt.model.ForecastU;
|
||||
import com.gunshi.project.xyt.model.ForecastUseparam;
|
||||
import com.gunshi.project.xyt.service.ForecastUService;
|
||||
import com.gunshi.project.xyt.service.ForecastUseparamService;
|
||||
import com.gunshi.project.xyt.validate.markers.Insert;
|
||||
import com.gunshi.project.xyt.validate.markers.Update;
|
||||
|
|
@ -23,9 +27,13 @@ import org.springframework.web.bind.annotation.RequestMapping;
|
|||
import org.springframework.web.bind.annotation.RestController;
|
||||
|
||||
import java.io.Serializable;
|
||||
import java.math.BigDecimal;
|
||||
import java.util.ArrayList;
|
||||
import java.util.Date;
|
||||
import java.util.List;
|
||||
import java.util.Map;
|
||||
import java.util.Objects;
|
||||
import java.util.stream.Collectors;
|
||||
|
||||
/**
|
||||
* 描述: 预报_通用参数管理
|
||||
|
|
@ -40,6 +48,9 @@ public class ForecastUseparamController {
|
|||
@Autowired
|
||||
private ForecastUseparamService service;
|
||||
|
||||
@Autowired
|
||||
private ForecastUService forecastUService;
|
||||
|
||||
|
||||
@Operation(summary = "新增")
|
||||
@PostMapping("/insert")
|
||||
|
|
@ -59,6 +70,37 @@ public class ForecastUseparamController {
|
|||
dto.setChtm(oldData.getChtm());
|
||||
dto.setUpdateTime(new Date());
|
||||
boolean result = service.updateById(dto);
|
||||
// 修改计算单位线的参数时，需要重新计算一次
|
||||
if (result && (dto.getParamCode().equals("cymj") || dto.getParamCode().equals("hdc") || dto.getParamCode().equals("hdpj") || dto.getParamCode().equals("dt") || dto.getParamCode().equals("h")) && !dto.getParamValue().equals(oldData.getParamValue())) {
|
||||
// 承雨面积areaF："cymj", 河道长lengthL："hdc", 河道坡降j："hdpj", 时段∆t："dt", h毫米净雨："h"
|
||||
QueryWrapper<ForecastUseparam> qw = new QueryWrapper<ForecastUseparam>().in("param_code", "cymj", "hdc", "hdpj", "dt", "h");
|
||||
List<Map<String, Object>> uMaps = service.listMaps(qw);
|
||||
Map<String, String> uMap = uMaps.stream().collect(Collectors.toMap(map -> (String) map.get("param_code"), map -> (String) map.get("param_value")));
|
||||
if (uMap.containsKey("cymj") && StringUtils.isNotEmpty(uMap.get("cymj"))
|
||||
&& uMap.containsKey("hdc") && StringUtils.isNotEmpty(uMap.get("hdc"))
|
||||
&& uMap.containsKey("hdpj") && StringUtils.isNotEmpty(uMap.get("hdpj"))
|
||||
&& uMap.containsKey("dt") && StringUtils.isNotEmpty(uMap.get("dt"))) {
|
||||
List<BigDecimal> uParam = service.calcU(new BigDecimal(uMap.get("cymj")),
|
||||
new BigDecimal(uMap.get("hdc")),
|
||||
new BigDecimal(uMap.get("hdpj")),
|
||||
Double.valueOf(uMap.get("dt")),
|
||||
StringUtils.isNotEmpty(uMap.get("h")) ? new BigDecimal(uMap.get("h")) : BigDecimal.ONE, StringUtils.isNotEmpty(uMap.get("swfq")) ? uMap.get("swfq") : "I");
|
||||
if (CollectionUtils.isNotEmpty(uParam)) {
|
||||
forecastUService.remove(new UpdateWrapper<>());
|
||||
List<ForecastU> uList = new ArrayList<>();
|
||||
Date date = new Date();
|
||||
for (BigDecimal u : uParam) {
|
||||
ForecastU forecastU = new ForecastU();
|
||||
forecastU.setId(IdWorker.getId());
|
||||
forecastU.setUValue(u);
|
||||
forecastU.setChtm(date);
|
||||
forecastU.setUpdateTime(date);
|
||||
uList.add(forecastU);
|
||||
}
|
||||
forecastUService.saveBatch(uList);
|
||||
}
|
||||
}
|
||||
}
|
||||
return R.ok(result ? dto : null);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -761,106 +761,4 @@ public class ForecastResultsService extends ServiceImpl<ForecastResultsMapper, F
|
|||
jsonBuilder.append("}"); // 闭合JSON对象
|
||||
return jsonBuilder.toString();
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* @description: 计算时间单位线u参数
|
||||
* @param areaF 流域面积 k㎡
|
||||
* @param lengthL 主河道长度 km
|
||||
* @param j 主河道加权平均比降 千分率‰
|
||||
* @param dt 时段∆t，时间间隔
|
||||
* @param tl 应为时间跨度，例如72小时
|
||||
* @return: java.util.List<java.math.BigDecimal>
|
||||
* @auther: cxw
|
||||
* @date: 2024-11-06, 周三, 10:16:26
|
||||
*/
|
||||
private List<BigDecimal> calcU(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, double dt, BigDecimal tl){
|
||||
// 计算m1、n
|
||||
Map<String, BigDecimal> m1AndNMap = calcM1AndN(areaF, lengthL, j);
|
||||
BigDecimal m1 = m1AndNMap.get("m1");
|
||||
BigDecimal n = m1AndNMap.get("n");
|
||||
// k = m1 / n。如果m1有非线性改正，则为改正后的值。算例保留1位小数
|
||||
BigDecimal k = m1.divide(n, 1, BigDecimal.ROUND_HALF_UP);
|
||||
/*
|
||||
* 瞬时单位线u(0,t)=(1/(k*Γ(n)))*((t/k)^(n-1))*e^(-t/k)
|
||||
*/
|
||||
// 计算伽马函数
|
||||
BigDecimal nGamma = BigDecimal.ONE.divide(k.multiply(BigdecimalUtil.bigDecimalGamma(n)), 4, BigDecimal.ROUND_HALF_UP);
|
||||
List<BigDecimal> uList = new ArrayList<>();
|
||||
for (int i = 1; i <= 11; i++) {
|
||||
BigDecimal t = BigDecimal.valueOf(dt * i);
|
||||
if (t.compareTo(tl) > 0) {
|
||||
break;
|
||||
}
|
||||
// 计算t、k、n
|
||||
BigDecimal tkn = BigdecimalUtil.bigDecimalExponentiation(t.divide(k, 4, BigDecimal.ROUND_HALF_UP), n.subtract(BigDecimal.ONE).doubleValue()).setScale(2, BigDecimal.ROUND_HALF_UP);
|
||||
// 计算 e 的指数
|
||||
BigDecimal eToThePower = BigdecimalUtil.bigDecimalExponentiation(BigDecimal.valueOf(Math.E), t.divide(k, 4, BigDecimal.ROUND_HALF_UP).multiply(BigDecimal.ONE.negate()).setScale(2, BigDecimal.ROUND_HALF_UP).doubleValue());
|
||||
// 计算u(0,t)
|
||||
BigDecimal u0t = nGamma.multiply(tkn).multiply(eToThePower).setScale(2, BigDecimal.ROUND_HALF_UP);
|
||||
uList.add(u0t);
|
||||
}
|
||||
System.out.println(uList);
|
||||
// TODO S(t)=0∫t*u(0,t)dt为参数n、k的函数，有S(t)查算表？
|
||||
|
||||
// TODO 时段单位线公式u(∆t,t)=[S(t)-S(t-∆t)]
|
||||
|
||||
// TODO h毫米净雨时段单位线q(∆t,t)=Fh/(3.6*∆t)*[u(∆t,t)]
|
||||
|
||||
|
||||
return new ArrayList<>();
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* @description: 计算m1和n参数
|
||||
* @param areaF
|
||||
* @param lengthL
|
||||
* @param j
|
||||
* @return: void
|
||||
* @auther: cxw
|
||||
* @date: 2024-11-06, 周三, 13:16:24
|
||||
*/
|
||||
private Map<String, BigDecimal> calcM1AndN(BigDecimal areaF, BigDecimal lengthL, BigDecimal j) {
|
||||
BigDecimal m1 = BigDecimal.ZERO;
|
||||
BigDecimal n = BigDecimal.ZERO;
|
||||
// 所属片区（全省山丘区瞬时单位线分三个片进行参数的地区综合，第一片包括水文分区 1、2、4区即京广线两侧及鄂东黄冈、咸宁地区一带,(江汉平原湖区在外);第二片包括水文分区6、8、9、11区即鄂北，鄂西北及宜昌地区长江以北一带;第三片包括7、10水文分区即清江流域，恩施地区）
|
||||
String slice = "I";
|
||||
if (slice.toUpperCase().equals("I")) {// I片（1、2、4区）
|
||||
if (areaF.compareTo(new BigDecimal("30")) > 0) {// areaF > 30
|
||||
// m1 = 0.82*F^0.29*L^0.23*j^-0.20
|
||||
m1 = new BigDecimal("0.82").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.29)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.23)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.20));
|
||||
} else if (areaF.compareTo(new BigDecimal("30")) <= 0) {// areaF <= 30
|
||||
// m1 = 1.38*F^0.27*L^0.216*j^-0.185
|
||||
m1 = new BigDecimal("1.38").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.27)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.216)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.185));
|
||||
}
|
||||
if (j.compareTo(new BigDecimal("5")) > 0) {// j > 5
|
||||
// n = 0.34*F^0.35*j^0.1
|
||||
n = new BigDecimal("0.34").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.35)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1));
|
||||
} else if (j.compareTo(new BigDecimal("5")) <= 0) {// j <= 5
|
||||
// n = 1.04*F^0.3*j^0.1
|
||||
n = new BigDecimal("1.04").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1));
|
||||
}
|
||||
} else if (slice.toUpperCase().equals("II")) {// II片（6、8、9、11区）
|
||||
// m1 = 1.64*F^0.231*L^0.131*j^-0.08
|
||||
m1 = new BigDecimal("1.64").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.231)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.131)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.08));
|
||||
// n = 0.529*F^0.25*j^0.20
|
||||
n = new BigDecimal("0.529").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.25)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.20));
|
||||
} else if (slice.toUpperCase().equals("III")) {// III片（7、10区）
|
||||
// m1 = 0.8*F^0.3*L^0.1*j^-0.06
|
||||
m1 = new BigDecimal("0.8").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.1)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.06));
|
||||
// n = 0.69*F^0.224*j^0.092
|
||||
n = new BigDecimal("0.69").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.224)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.092));
|
||||
}
|
||||
// TODO 特殊流域适用的公式
|
||||
|
||||
// TODO 当设计雨强超过10mm/h时，除熔岩地区外一般应作参数m1外延的非线性改正
|
||||
|
||||
BigDecimal finalM1 = m1.setScale(1, BigDecimal.ROUND_HALF_UP);
|
||||
BigDecimal finalN = n.setScale(1, BigDecimal.ROUND_HALF_UP);
|
||||
return new HashMap<>() {{
|
||||
put("m1", finalM1);
|
||||
put("n", finalN);
|
||||
}};
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -3,10 +3,18 @@ package com.gunshi.project.xyt.service;
|
|||
import com.baomidou.mybatisplus.extension.service.impl.ServiceImpl;
|
||||
import com.gunshi.project.xyt.mapper.ForecastUseparamMapper;
|
||||
import com.gunshi.project.xyt.model.ForecastUseparam;
|
||||
import com.gunshi.project.xyt.util.BigdecimalUtil;
|
||||
import lombok.extern.slf4j.Slf4j;
|
||||
import org.springframework.stereotype.Service;
|
||||
import org.springframework.transaction.annotation.Transactional;
|
||||
|
||||
import java.math.BigDecimal;
|
||||
import java.util.ArrayList;
|
||||
import java.util.HashMap;
|
||||
import java.util.List;
|
||||
import java.util.Map;
|
||||
import java.util.stream.Collectors;
|
||||
|
||||
/**
|
||||
* 描述: 预报_通用参数管理
|
||||
* author: cxw
|
||||
|
|
@ -18,6 +26,122 @@ import org.springframework.transaction.annotation.Transactional;
|
|||
public class ForecastUseparamService extends ServiceImpl<ForecastUseparamMapper, ForecastUseparam>
|
||||
{
|
||||
|
||||
/**
|
||||
* @description: 计算时间单位线u参数
|
||||
* @param areaF 流域面积 k㎡（承雨面积）
|
||||
* @param lengthL 主河道长度 km
|
||||
* @param j 主河道加权平均比降 千分率‰（河道坡降）
|
||||
* @param dt 时段∆t，时间间隔，例如0.5
|
||||
* @param h h毫米净雨，例如1、10等，默认为1
|
||||
* @param slice 水文分区
|
||||
* @return: java.util.List<java.math.BigDecimal>
|
||||
* @auther: cxw
|
||||
* @date: 2024-11-06, 周三, 10:16:26
|
||||
*/
|
||||
public List<BigDecimal> calcU(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, double dt, BigDecimal h, String slice) {
|
||||
// 计算m1、n
|
||||
Map<String, BigDecimal> m1AndNMap = calcM1AndN(areaF, lengthL, j, slice);
|
||||
BigDecimal m1 = m1AndNMap.get("m1");
|
||||
BigDecimal n = m1AndNMap.get("n");
|
||||
// k = m1 / n。如果m1有非线性改正，则为改正后的值。算例保留1位小数
|
||||
BigDecimal k = m1.divide(n, 3, BigDecimal.ROUND_HALF_UP);
|
||||
/*
|
||||
* 瞬时单位线u(0,t)=(1/(k*Γ(n)))*((t/k)^(n-1))*e^(-t/k)
|
||||
*/
|
||||
List<Map<String, BigDecimal>> mapList = new ArrayList<>();
|
||||
Map<String, BigDecimal> initMap = new HashMap<>();
|
||||
initMap.put("t", BigDecimal.ZERO);
|
||||
initMap.put("tDevideK", BigDecimal.ZERO);
|
||||
initMap.put("st", BigDecimal.ZERO);
|
||||
initMap.put("utt", BigDecimal.ZERO);
|
||||
initMap.put("qtt", BigDecimal.ZERO);
|
||||
mapList.add(initMap);
|
||||
for (int i = 1; i <= 11; i++) {
|
||||
Map<String, BigDecimal> result = new HashMap<>();
|
||||
// 表1-5地表径流过程计算表：t
|
||||
BigDecimal t = BigDecimal.valueOf(dt * i);
|
||||
result.put("t", t);
|
||||
|
||||
// 表1-5地表径流过程计算表：t/k
|
||||
BigDecimal x = t.divide(k, 10, BigDecimal.ROUND_HALF_UP);// x
|
||||
result.put("tDevideK", x);
|
||||
|
||||
BigDecimal alpha = n;// α
|
||||
BigDecimal beta = BigDecimal.ONE;// β
|
||||
// 表1-5地表径流过程计算表：s(t)
|
||||
BigDecimal gammaDist = BigDecimal.valueOf(BigdecimalUtil.gammaDist(x.doubleValue(), alpha.doubleValue(), beta.doubleValue(), true)).setScale(10, BigDecimal.ROUND_HALF_UP);// 伽马分布
|
||||
result.put("st", gammaDist);
|
||||
// 表1-5地表径流过程计算表：s(t-∆t)。stt为上一行的st
|
||||
BigDecimal stt = BigDecimal.ZERO;
|
||||
stt = mapList.get(i - 1).get("st").setScale(10, BigDecimal.ROUND_HALF_UP);
|
||||
result.put("stt", stt);
|
||||
|
||||
// 表1-5地表径流过程计算表：u(∆t,t)。utt=st-stt
|
||||
BigDecimal utt = gammaDist.subtract(stt).setScale(10, BigDecimal.ROUND_HALF_UP);
|
||||
result.put("utt", utt);
|
||||
|
||||
// 表1-5地表径流过程计算表：q(∆t,t)。qtt=Fh*utt/(3.6*∆t)
|
||||
BigDecimal qtt = areaF.multiply(h).multiply(utt).divide(new BigDecimal(3.6).multiply(t), 3, BigDecimal.ROUND_HALF_UP);
|
||||
result.put("qtt", qtt);
|
||||
mapList.add(result);
|
||||
}
|
||||
mapList.remove(0);// 去除initMap
|
||||
List<BigDecimal> uList = mapList.stream().map(map -> map.get("qtt")).collect(Collectors.toList());
|
||||
return uList;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* @description: 计算m1和n参数
|
||||
* @param areaF
|
||||
* @param lengthL
|
||||
* @param j
|
||||
* @param slice
|
||||
* @return: void
|
||||
* @auther: cxw
|
||||
* @date: 2024-11-06, 周三, 13:16:24
|
||||
*/
|
||||
private Map<String, BigDecimal> calcM1AndN(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, String slice) {
|
||||
BigDecimal m1 = BigDecimal.ZERO;
|
||||
BigDecimal n = BigDecimal.ZERO;
|
||||
// 所属片区（全省山丘区瞬时单位线分三个片进行参数的地区综合，第一片包括水文分区 1、2、4区即京广线两侧及鄂东黄冈、咸宁地区一带,(江汉平原湖区在外);第二片包括水文分区6、8、9、11区即鄂北，鄂西北及宜昌地区长江以北一带;第三片包括7、10水文分区即清江流域，恩施地区）
|
||||
if (slice.toUpperCase().equals("I")) {// I片（1、2、4区）
|
||||
if (areaF.compareTo(new BigDecimal("30")) > 0) {// areaF > 30
|
||||
// m1 = 0.82*F^0.29*L^0.23*j^-0.20
|
||||
m1 = new BigDecimal("0.82").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.29)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.23)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.20));
|
||||
} else if (areaF.compareTo(new BigDecimal("30")) <= 0) {// areaF <= 30
|
||||
// m1 = 1.38*F^0.27*L^0.216*j^-0.185
|
||||
m1 = new BigDecimal("1.38").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.27)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.216)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.185));
|
||||
}
|
||||
if (j.compareTo(new BigDecimal("5")) > 0) {// j > 5
|
||||
// n = 0.34*F^0.35*j^0.1
|
||||
n = new BigDecimal("0.34").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.35)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1));
|
||||
} else if (j.compareTo(new BigDecimal("5")) <= 0) {// j <= 5
|
||||
// n = 1.04*F^0.3*j^0.1
|
||||
n = new BigDecimal("1.04").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1));
|
||||
}
|
||||
} else if (slice.toUpperCase().equals("II")) {// II片（6、8、9、11区）
|
||||
// m1 = 1.64*F^0.231*L^0.131*j^-0.08
|
||||
m1 = new BigDecimal("1.64").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.231)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.131)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.08));
|
||||
// n = 0.529*F^0.25*j^0.20
|
||||
n = new BigDecimal("0.529").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.25)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.20));
|
||||
} else if (slice.toUpperCase().equals("III")) {// III片（7、10区）
|
||||
// m1 = 0.8*F^0.3*L^0.1*j^-0.06
|
||||
m1 = new BigDecimal("0.8").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.1)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.06));
|
||||
// n = 0.69*F^0.224*j^0.092
|
||||
n = new BigDecimal("0.69").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.224)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.092));
|
||||
}
|
||||
// TODO 特殊流域适用的公式
|
||||
|
||||
// TODO 当设计雨强超过10mm/h时，除熔岩地区外一般应作参数m1外延的非线性改正
|
||||
|
||||
BigDecimal finalM1 = m1.setScale(3, BigDecimal.ROUND_HALF_UP);
|
||||
BigDecimal finalN = n.setScale(3, BigDecimal.ROUND_HALF_UP);
|
||||
return new HashMap<>() {{
|
||||
put("m1", finalM1);
|
||||
put("n", finalN);
|
||||
}};
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,7 @@
|
|||
package com.gunshi.project.xyt.util;
|
||||
|
||||
import org.apache.commons.math3.distribution.GammaDistribution;
|
||||
|
||||
import java.math.BigDecimal;
|
||||
import java.math.MathContext;
|
||||
import java.math.RoundingMode;
|
||||
|
|
@ -34,18 +36,12 @@ public class BigdecimalUtil {
|
|||
* @auther: cxw
|
||||
* @date: 2024-11-06, 周三, 11:31:24
|
||||
*/
|
||||
public static BigDecimal bigDecimalGamma(BigDecimal x) {
|
||||
BigDecimal result = BigDecimal.valueOf(gamma(x.doubleValue()));
|
||||
return result.setScale(x.scale(), RoundingMode.HALF_UP);
|
||||
}
|
||||
public static double gamma(double x) {
|
||||
return Math.exp(logGamma(x));
|
||||
}
|
||||
public static double logGamma(double x) {
|
||||
double tmp = (x - 0.5) * Math.log(x + 4.5) - (x + 4.5);
|
||||
double ser = 1.0 + 76.18009173 / (x + 0) - 86.50532033 / (x + 1)
|
||||
+ 24.01409822 / (x + 2) - 1.231739516 / (x + 3)
|
||||
+ 0.00120858003 / (x + 4) - 0.00000536382 / (x + 5);
|
||||
return tmp + Math.log(ser * Math.sqrt(2 * Math.PI));
|
||||
public static double gammaDist(double x, double alpha, double beta, boolean cumulative) {
|
||||
GammaDistribution gammaDist = new GammaDistribution(alpha, 1/beta);
|
||||
if (cumulative) {
|
||||
return gammaDist.cumulativeProbability(x);
|
||||
} else {
|
||||
return gammaDist.density(x);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 200 KiB |
Binary file not shown.
Loading…
Reference in New Issue