risk-assessment

name	risk-assessment
version	1.0.0
category	engineering
description	Risk Assessment Skill

Risk Assessment Skill

name: risk-assessment
version: 1.0.0
category: sme
tags: [risk, probabilistic, monte-carlo, reliability, uncertainty, sensitivity-analysis, marine-safety]
created: 2026-01-06
updated: 2026-01-06
author: Claude
description: |
  Expert risk assessment and probabilistic analysis for marine and offshore
  operations. Includes Monte Carlo simulations, reliability calculations,
  uncertainty quantification, and decision-making under risk.

When to Use This Skill

Use this skill when you need to:

• Perform Monte Carlo simulations for uncertainty quantification
• Calculate system reliability and failure probabilities
• Conduct sensitivity analysis to identify critical parameters
• Create risk matrices for hazard assessment
• Perform probabilistic design and analysis
• Quantify uncertainties in marine operations
• Make decisions under uncertainty with risk metrics
• Validate designs against reliability targets

Core Knowledge Areas

1. Monte Carlo Simulation

Basic Monte Carlo framework:

import numpy as np
from scipy import stats
from dataclasses import dataclass
from typing import Callable, Dict, List, Tuple, Optional
import pandas as pd

@dataclass
class RandomVariable:
    """Statistical distribution for a random variable."""
    name: str
    distribution: str  # 'normal', 'lognormal', 'uniform', 'weibull', etc.
    parameters: dict  # Distribution-specific parameters

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples from distribution.

        Args:
            size: Number of samples

        Returns:
            Array of random samples

        Example:
            >>> rv = RandomVariable(
            ...     name='wave_height',
            ...     distribution='weibull',
            ...     parameters={'c': 2.0, 'scale': 3.5}
            ... )
            >>> samples = rv.sample(1000)
        """
        if self.distribution == 'normal':
            return np.random.normal(
                loc=self.parameters['mean'],
                scale=self.parameters['std'],
                size=size
            )
        elif self.distribution == 'lognormal':
            return np.random.lognormal(
                mean=.parameters[],
                sigma=.parameters[],
                size=size
            )
         .distribution == :
             np.random.uniform(
                low=.parameters[],
                high=.parameters[],
                size=size
            )
         .distribution == :
             stats.weibull_min.rvs(
                c=.parameters[],
                scale=.parameters[],
                size=size
            )
         .distribution == :
             np.random.exponential(
                scale=.parameters[],
                size=size
            )
        :
             ValueError()

 () -> [, np.ndarray]:
    
     seed   :
        np.random.seed(seed)

    
    samples = {}
     rv  random_variables:
        samples[rv.name] = rv.sample(n_samples)

    
    outputs = model(samples)

    
    results = {**samples, : outputs}

     results

 () -> :
    
    stats_dict = {
        : np.mean(data),
        : np.std(data),
        : np.(data),
        : np.(data),
        : np.median(data),
        : np.std(data) / np.mean(data)  
    }

    
     p  percentiles:
        stats_dict[] = np.percentile(data, p)

     stats_dict

2. Reliability Analysis

Calculate reliability and failure probability:

def calculate_reliability(
    response_data: np.ndarray,
    limit_state: float,
    mode: str = 'less_than'
) -> dict:
    """
    Calculate reliability from Monte Carlo results.

    Args:
        response_data: Array of response values
        limit_state: Limit state threshold
        mode: 'less_than' or 'greater_than'

    Returns:
        Dictionary with reliability metrics

    Example:
        >>> # Mooring tension should be < 8000 kN
        >>> reliability = calculate_reliability(
        ...     response_data=results['output'],
        ...     limit_state=8000,
        ...     mode='less_than'
        ... )
        >>> print(f"Reliability: {reliability['reliability']:.4f}")
        >>> print(f"Probability of failure: {reliability['pf']:.6f}")
    """
    n_samples = len(response_data)

    if mode == 'less_than':
        # Failure when response >= limit_state
        failures = response_data >= limit_state
    elif mode == 'greater_than':
        # Failure when response <= limit_state
        failures = response_data <= limit_state
    else:
        raise ValueError(f"Unknown mode: {mode}")

    n_failures = np.sum(failures)
    pf = n_failures / n_samples  # Probability of failure
    reliability = 1 - pf

    # Reliability index (beta)
    # β = -Φ^(-1)(Pf) where Φ is standard normal CDF
    if pf > 0 and pf < 1:
        beta = -stats.norm.ppf(pf)
    elif pf == 0:
        beta = np.inf
    else:  # pf == 1
        beta = -np.inf

    return {
        'n_samples': n_samples,
        'n_failures': n_failures,
        : pf,
        : reliability,
        : beta,
        : reliability >=   
    }

 () -> :
    
    
    R_sys = np.prod(component_reliabilities)
     R_sys

 () -> :
    
    
    pf_components = [ - R  R  component_reliabilities]
    pf_sys = np.prod(pf_components)
    R_sys =  - pf_sys
     R_sys

 () -> :
    
    
    

    numerator = np.dot(gradient, mean)
    denominator = np.sqrt(
        np.dot(gradient, np.dot(cov_matrix, gradient))
    )

    beta = numerator / denominator

    
    pf = stats.norm.cdf(-beta)

     {
        : beta,
        : pf,
        :  - pf
    }

3. Sensitivity Analysis

Identify critical parameters:

def sensitivity_analysis_correlation(
    inputs: Dict[str, np.ndarray],
    output: np.ndarray
) -> pd.DataFrame:
    """
    Sensitivity analysis using correlation coefficients.

    Args:
        inputs: Dictionary of input variable samples
        output: Output variable samples

    Returns:
        DataFrame with correlation coefficients sorted by absolute value

    Example:
        >>> sensitivity = sensitivity_analysis_correlation(
        ...     inputs={
        ...         'wave_height': results['wave_height'],
        ...         'current_speed': results['current_speed'],
        ...         'wind_speed': results['wind_speed']
        ...     },
        ...     output=results['output']
        ... )
        >>> print(sensitivity)
    """
    correlations = {}

    for var_name, var_samples in inputs.items():
        # Pearson correlation coefficient
        corr = np.corrcoef(var_samples, output)[0, 1]
        correlations[var_name] = corr

    # Create DataFrame
    df = pd.DataFrame.from_dict(
        correlations,
        orient='index',
        columns=['correlation']
    )

    # Add absolute value and rank
    df['abs_correlation'] = df['correlation'].abs()
    df = df.sort_values('abs_correlation', ascending=False)
    df['rank'] = range(1, len(df) + 1)

    return df

def sensitivity_analysis_variance(
    inputs: Dict[str, np.ndarray],
    output: np.ndarray,
    n_partitions: int = 10
) -> pd.DataFrame:
    """
    Variance-based sensitivity analysis (Sobol indices approximation).

    Args:
        inputs: Dictionary of input variable samples
        output: Output variable samples
        n_partitions: Number of partitions for conditional variance

    Returns:
        DataFrame with sensitivity indices

    Example:
        >>> sensitivity = sensitivity_analysis_variance(
        ...     inputs=results,
        ...     output=results['output']
        ... )
    """
    total_variance = np.var(output)

    sensitivity_indices = {}

     var_name, var_samples  inputs.items():
         var_name == :
            

        
        percentiles = np.linspace(, , n_partitions + )
        bins = np.percentile(var_samples, percentiles)

        
        conditional_means = []
         i  (n_partitions):
            mask = (var_samples >= bins[i]) & (var_samples < bins[i + ])
             np.(mask) > :
                conditional_means.append(np.mean(output[mask]))

        
        variance_of_means = np.var(conditional_means)

        
        S1 = variance_of_means / total_variance  total_variance >   

        sensitivity_indices[var_name] = {
            : S1
        }

    
    df = pd.DataFrame.from_dict(sensitivity_indices, orient=)
    df = df.sort_values(, ascending=)
    df[] = (, (df) + )

     df

 () -> pd.DataFrame:
    
    results = []

     var_name, var_values  variables.items():
        
        

        inputs_low = {k: v[]  k, v  variables.items()}
        inputs_low[var_name] = var_values[]
        output_low = model(inputs_low)

        inputs_high = {k: v[]  k, v  variables.items()}
        inputs_high[var_name] = var_values[]
        output_high = model(inputs_high)

        
        swing = (output_high - output_low)

        results.append({
            : var_name,
            : output_low,
            : output_high,
            : swing,
            : var_values[],
            : var_values[]
        })

    df = pd.DataFrame(results)
    df = df.sort_values(, ascending=)
    df[] = (, (df) + )

     df

4. Risk Matrices and Hazard Assessment

from enum import Enum

class Severity(Enum):
    """Consequence severity levels."""
    NEGLIGIBLE = 1
    MINOR = 2
    MODERATE = 3
    MAJOR = 4
    CATASTROPHIC = 5

class Likelihood(Enum):
    """Event likelihood levels."""
    RARE = 1
    UNLIKELY = 2
    POSSIBLE = 3
    LIKELY = 4
    ALMOST_CERTAIN = 5

class RiskLevel(Enum):
    """Risk level categories."""
    LOW = 1
    MEDIUM = 2
    HIGH = 3
    VERY_HIGH = 4

@dataclass
class Hazard:
    """Hazard definition."""
    id: str
    description: str
    severity: Severity
    likelihood: Likelihood
    existing_controls: List[str]
    residual_risk: Optional[RiskLevel] = None

def calculate_risk_level(
    severity: Severity,
    likelihood: Likelihood
) -> RiskLevel:
    """
    Calculate risk level from severity and likelihood.

    Uses 5x5 risk matrix.

    Args:
        severity: Consequence severity
        likelihood: Event likelihood

    Returns:
        Risk level

    Example:
        >>> risk = calculate_risk_level(
        ...     Severity.MAJOR,
        ...     Likelihood.POSSIBLE
        ... )
        >>> print(risk)  # RiskLevel.HIGH
    """
    
    
    
    risk_matrix = np.array([
        [, , , , ],  
        [, , , , ],  
        [, , , , ],  
        [, , , , ],  
        [, , , , ]   
    ])

    risk_value = risk_matrix[likelihood.value - , severity.value - ]

    
     risk_value == :
         RiskLevel.LOW
     risk_value == :
         RiskLevel.MEDIUM
     risk_value == :
         RiskLevel.HIGH
    :  
         RiskLevel.VERY_HIGH

 () -> pd.DataFrame:
    
    results = []

     hazard  hazards:
        risk_level = calculate_risk_level(hazard.severity, hazard.likelihood)

        results.append({
            : hazard.,
            : hazard.description,
            : hazard.severity.name,
            : hazard.likelihood.name,
            : risk_level.name,
            : .join(hazard.existing_controls)
        })

    df = pd.DataFrame(results)

    
    risk_order = {
        : ,
        : ,
        : ,
        : 
    }
    df[] = df[].(risk_order)
    df = df.sort_values(, ascending=)
    df = df.drop(, axis=)

     df

5. Extreme Value Analysis

from scipy.stats import genextreme

def fit_extreme_value_distribution(
    data: np.ndarray,
    method: str = 'gev'
) -> dict:
    """
    Fit extreme value distribution to data.

    Args:
        data: Array of extreme values (e.g., annual maxima)
        method: 'gev' (Generalized Extreme Value) or 'gumbel'

    Returns:
        Dictionary with fitted parameters and statistics

    Example:
        >>> # Annual maximum wave heights
        >>> annual_max_Hs = np.array([8.5, 9.2, 7.8, 10.1, 8.9, ...])
        >>> fit = fit_extreme_value_distribution(annual_max_Hs, 'gev')
        >>> print(f"Shape: {fit['shape']:.3f}")
        >>> print(f"Location: {fit['location']:.2f}")
        >>> print(f"Scale: {fit['scale']:.2f}")
    """
    if method == 'gev':
        # Fit GEV distribution
        shape, location, scale = genextreme.fit(data)

        # Calculate return values
        return_periods = [10, 50, 100, 1000, 10000]  # years
        return_values = {}

        for T in return_periods:
            # Return level for return period T
            p = 1 - 1/T
            x_T = genextreme.ppf(p, shape, loc=location, scale=scale)
            return_values[f'{T}yr'] = x_T

        return {
            'shape': shape,
            'location': location,
            'scale': scale,
            'return_values': return_values
        }

    elif method == 'gumbel':
        # Gumbel is special case of GEV with shape=0
        location, scale = stats.gumbel_r.fit(data)

        return_periods = [10, , , , ]
        return_values = {}

         T  return_periods:
            p =  - /T
            x_T = stats.gumbel_r.ppf(p, loc=location, scale=scale)
            return_values[] = x_T

         {
            : ,  
            : location,
            : scale,
            : return_values
        }

    :
         ValueError()

Complete Examples

Example 1: Complete Mooring System Risk Assessment

import numpy as np
from pathlib import Path

def complete_mooring_risk_assessment(
    design_parameters: dict,
    environmental_parameters: dict,
    n_simulations: int = 10000
) -> dict:
    """
    Complete probabilistic risk assessment for mooring system.

    Example:
        >>> design = {
        ...     'n_lines': 8,
        ...     'line_capacity': 10000,  # kN
        ...     'safety_factor': 2.5
        ... }
        >>> environment = {
        ...     'wave_height': {'distribution': 'weibull', 'c': 2.0, 'scale': 3.5},
        ...     'wave_period': {'distribution': 'normal', 'mean': 10.0, 'std': 2.0},
        ...     'current_speed': {'distribution': 'normal', 'mean': 1.0, 'std': 0.3}
        ... }
        >>> results = complete_mooring_risk_assessment(
        ...     design,
        ...     environment,
        ...     n_simulations=10000
        ... )
    """
    print("="*70)
    print("MOORING SYSTEM RISK ASSESSMENT")
    print("="*70)

    # Step 1: Define random variables
    print("\n[Step 1/6] Defining random variables...")

    random_vars = []
    for var_name, var_params in environment_parameters.items():
        rv = RandomVariable(
            name=var_name,
            distribution=var_params['distribution'],
            parameters={k: v for k, v in var_params.items() if k != 'distribution'}
        )
        random_vars.append(rv)
        print(f"  {var_name}: {var_params['distribution']}")

    
    ()

     ():
        
        Hs = samples[]
        Tp = samples[]
        V_c = samples[]

        
        tension = (
             +  
             * Hs +  
             * Hs** / Tp +  
             * V_c  
        )

         tension

    
    ()

    mc_results = monte_carlo_simulation(
        model=mooring_tension_model,
        random_variables=random_vars,
        n_samples=n_simulations,
        seed=
    )

    
    stats = calculate_statistics(mc_results[])

    ()
    ()
    ()
    ()

    
    ()

    
    design_limit = design_parameters[] / design_parameters[]

    reliability_result = calculate_reliability(
        response_data=mc_results[],
        limit_state=design_limit,
        mode=
    )

    ()
    ()
    ()
    ()

    
    ()

    sensitivity = sensitivity_analysis_correlation(
        inputs={k: v  k, v  mc_results.items()  k != },
        output=mc_results[]
    )

    ()
     idx, row  sensitivity.iterrows():
        ()

    
    ()

    
    component_reliability = reliability_result[]

    
    R_all_intact = system_reliability_series(
        [component_reliability] * design_parameters[]
    )

    
    R_total_failure =  - ( - component_reliability)**design_parameters[]

    ()
    ()
    ()

    
    summary = {
        : stats,
        : reliability_result,
        : sensitivity,
        : {
            : R_all_intact,
            : R_total_failure
        }
    }

    ( + *)
    ()
    (*)

     summary


design_parameters = {
    : ,
    : ,
    : 
}

environmental_parameters = {
    : {: , : , : },
    : {: , : , : },
    : {: , : , : }
}

risk_assessment = complete_mooring_risk_assessment(
    design_parameters,
    environmental_parameters,
    n_simulations=
)

Best Practices

1. Sample Size Selection

def determine_sample_size(
    target_pf: float,
    confidence_level: float = 0.95
) -> int:
    """
    Determine required sample size for target failure probability.

    Rule of thumb: N ≈ 10/Pf for reasonable confidence

    Args:
        target_pf: Target probability of failure
        confidence_level: Confidence level

    Returns:
        Recommended sample size

    Example:
        >>> n = determine_sample_size(target_pf=1e-4, confidence_level=0.95)
        >>> print(f"Recommended samples: {n}")
    """
    # Basic rule: need at least 10 failures
    # So N * Pf ≥ 10
    # N ≥ 10 / Pf

    n_basic = int(10 / target_pf)

    # For higher confidence, increase further
    if confidence_level >= 0.99:
        n_recommended = n_basic * 5
    elif confidence_level >= 0.95:
        n_recommended = n_basic * 3
    else:
        n_recommended = n_basic * 2

    return max(n_recommended, 1000)  # Minimum 1000 samples

2. Convergence Checking

def check_monte_carlo_convergence(
    data: np.ndarray,
    window_size: int = 1000
) -> dict:
    """
    Check if Monte Carlo simulation has converged.

    Args:
        data: Simulation output data
        window_size: Window size for moving average

    Returns:
        Dictionary with convergence metrics
    """
    n = len(data)

    # Calculate cumulative mean
    cumulative_mean = np.cumsum(data) / np.arange(1, n + 1)

    # Calculate moving coefficient of variation
    if n > window_size:
        moving_std = np.std(data[-window_size:])
        moving_mean = np.mean(data[-window_size:])
        cov = moving_std / moving_mean if moving_mean != 0 else 0
    else:
        cov = np.std(data) / np.mean(data) if np.mean(data) != 0 else 0

    # Convergence criterion: COV < 5%
    converged = cov < 0.05

    return {
        'converged': converged,
        'cov': cov,
        'final_mean': cumulative_mean[-1],
        'samples_used': n
    }

Resources

Textbooks

• Ang, A.H-S., Tang, W.H. (2007). Probability Concepts in Engineering
• Melchers, R.E., Beck, A.T. (2018). Structural Reliability Analysis and Prediction
• DNV (2021). DNVGL-RP-C205: Environmental Conditions and Environmental Loads

Standards

• DNV-RP-C205: Environmental conditions and environmental loads
• DNV-RP-C206: Fatigue methodology of offshore ships
• ISO 2394: General principles on reliability for structures
• API RP 2A: Planning, designing and constructing fixed offshore platforms

Software

• @RISK: Monte Carlo simulation add-in for Excel
• Crystal Ball: Oracle's risk analysis software
• OpenTURNS: Open source library for uncertainty quantification
• scipy.stats: Python statistical distributions

---

Use this skill for: Expert probabilistic risk assessment and reliability analysis for marine and offshore systems with comprehensive uncertainty quantification.